diff --git a/examples/2d/base/synthetic_solve.py b/examples/2d/base/synthetic_solve.py
new file mode 100644
index 0000000..6e58ce7
--- /dev/null
+++ b/examples/2d/base/synthetic_solve.py
@@ -0,0 +1,26 @@
+import numpy as np
+import fenics as fen
+from rrompy.hfengines.linear_problem.bidimensional import \
+ SyntheticBivariateEngine as SBE
+
+verb = 100
+
+kappa = 15. ** .5
+theta = - np.pi / 6.
+n = 30
+L = np.pi
+
+mu0 = [15. ** .5, 20. ** .5]
+mutar = [15. ** .5, 20. ** .5]
+
+solver = SBE(kappa = kappa, theta = theta, n = n, L = L,
+ mu0 = mu0, verbosity = verb)
+uh = solver.solve(mutar)[0]
+solver.plotmesh(figsize = (7.5, 4.5))
+fen.plot(fen.project(solver._above,
+ fen.FunctionSpace(solver.V.mesh(), "DG", 0)))
+print(solver.norm(uh))
+solver.plot(uh)
+#solver.outParaviewTimeDomain(uh, mutar[0], filename = 'out', folder = True,
+# forceNewFile = False)
+
diff --git a/examples/2d/pod/cookie_single_pod.py b/examples/2d/pod/cookie_single_pod.py
index 752762f..eac23b1 100644
--- a/examples/2d/pod/cookie_single_pod.py
+++ b/examples/2d/pod/cookie_single_pod.py
@@ -1,138 +1,137 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
CookieEngineSingle as CES
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
MN = 15
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if size == 1: # below
mu0 = [20 ** .5, 1. ** .5]
mutar = [20.5 ** .5, 1.05 ** .5]
murange = [[18.5 ** .5, .85 ** .5], [21.5 ** .5, 1.15 ** .5]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
kappa = 20. ** .5
theta = - np.pi / 6.
n = 40
Rad = 1.
L = np.pi
nX = 2
nY = 1
solver = CES(kappa = kappa, theta = theta, n = n, R = Rad, L = L, nX = nX,
nY = nY, mu0 = mu0, verbosity = verb)
rescaling = [lambda x: np.power(x, 2.)] * 2
rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', homogeneized = False,
what = "REAL")
approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL")
approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL")
# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- muInfVals, normEx, normApp, normErr = plotInfSet2(murange, murangeEff,
- approx, mu0, 40,
- [2., 2.], clip = clip)
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 40, [2., 2.], clip = clip)
diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index 64f053f..9c02db0 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,181 +1,191 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 7
+size = 3
show_sample = False
show_norm = True
-ignore_forcing = True
-ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
-MN = 20
+Delta = 0
+MN = 10
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
+PODTol = 1e-10
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
-#sampling = "random"
+sampling = "random"
+
+if samples == "distributed":
+ radial = 0
+# radial = "gaussian"
+# radial = "thinplate"
+# radial = "multiquadric"
+ rW0 = 10.
+ radialWeight = [r * rW0 for r in [40., .3]]
if size == 1: # below
mu0 = [40 ** .5, .4]
mutar = [45 ** .5, .4]
murange = [[30 ** .5, .3], [50 ** .5, .5]]
elif size == 2: # top
mu0 = [40 ** .5, .6]
mutar = [45 ** .5, .6]
murange = [[30 ** .5, .5], [50 ** .5, .7]]
elif size == 3: # interesting
mu0 = [40 ** .5, .5]
mutar = [45 ** .5, .5]
murange = [[30 ** .5, .3], [50 ** .5, .7]]
elif size == 4: # wide_low
mu0 = [40 ** .5, .2]
mutar = [45 ** .5, .2]
murange = [[10 ** .5, .1], [70 ** .5, .3]]
elif size == 5: # wide_hi
mu0 = [40 ** .5, .8]
mutar = [45 ** .5, .8]
murange = [[10 ** .5, .7], [70 ** .5, .9]]
elif size == 6: # top_zoom
mu0 = [50 ** .5, .8]
mutar = [55 ** .5, .8]
murange = [[40 ** .5, .7], [60 ** .5, .9]]
elif size == 7: # huge
mu0 = [50 ** .5, .5]
mutar = [55 ** .5, .5]
murange = [[10 ** .5, .2], [90 ** .5, .8]]
elif size == 100: # tiny
mu0 = [32.5 ** .5, .5]
mutar = [34 ** .5, .5]
murange = [[30 ** .5, .3], [35 ** .5, .7]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
H = 1.
L = .75
delta = .05
-n = 25
+n = 20
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescaling = [lambda x: np.power(x, 2.), lambda x: x]
rescalingInv = [lambda x: np.power(x, .5), lambda x: x]
if algo == "rational":
- params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN - Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
+ params['radialBasis'] = radial
+ params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
- params = {'R':R, 'S':S, 'POD':True}
+ params = {'R':(MN + 2 - Delta) * (MN + 1 - Delta) // 2, 'S':S,
+ 'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb, homogeneized = homogeneize)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = mutar[1]
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * L]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
approx.plotApprox(mutar, [warp, warpI], name = 'u_app',
homogeneized = False, what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF',
homogeneized = False, what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err',
homogeneized = False, what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res',
# homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar, homogeneized = homogeneize)
solNorm = approx.normHF(mutar, homogeneized = homogeneize)
resNorm = approx.normRes(mutar, homogeneized = homogeneize)
RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- muInfVals, normEx, normApp, normErr = plotInfSet2(murange, murangeEff,
- approx, mu0, 50,
- [2., 1.], clip = clip)
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 25, [2., 1.], clip = clip)
diff --git a/examples/2d/pod/fracture_pod_nodomain.py b/examples/2d/pod/fracture_pod_nodomain.py
index 428aaa4..7764988 100644
--- a/examples/2d/pod/fracture_pod_nodomain.py
+++ b/examples/2d/pod/fracture_pod_nodomain.py
@@ -1,176 +1,187 @@
import numpy as np
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 4
+size = 3
show_sample = True
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
-MN = 8
-R = MN + 1
+M = 5
+N = 6
+R = max(M, N) + 1
S = R
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
+if samples == "distributed":
+ radial = 0
+# radial = "gaussian"
+# radial = "thinplate"
+# radial = "multiquadric"
+ radialWeight = np.array([2.])
+
if size == 1: # below
mu0Aug = [40 ** .5, .4]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 2: # top
mu0Aug = [40 ** .5, .6]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 3: # interesting
mu0Aug = [40 ** .5, .5]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 4: # wide_low
mu0Aug = [40 ** .5, .2]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 5: # wide_hi
mu0Aug = [40 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 6: # top_zoom
mu0Aug = [50 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 55 ** .5
murange = [[40 ** .5], [60 ** .5]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5]]
H = 1.
L = .75
delta = .05
-n = 50
+n = 20
solver = MFEND(mu0 = mu0Aug, H = H, L = L, delta = delta, n = n,
verbosity = verb)
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
if algo == "rational":
- params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
+ params = {'N':N, 'M':M, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
- params['E'] = MN
+ params['E'] = max(M, N)
+ params['radialBasis'] = radial
+ params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
- params['S'] = [max(j, MN + 1) for j in params['S']]
+ params['S'] = [max(j, M + 1, N + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb, homogeneized = homogeneize)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = solver.lFrac
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * L]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
approx.plotApprox(mutar, [warp, warpI], name = 'u_app',
homogeneized = False, what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF',
homogeneized = False, what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err',
homogeneized = False, what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res',
# homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar, homogeneized = homogeneize)
solNorm = approx.normHF(mutar, homogeneized = homogeneize)
resNorm = approx.normRes(mutar, homogeneized = homogeneize)
RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet1
muZeroVals, Qvals = plotZeroSet1(murange, murangeEff, approx, mu0,
1000, 2.)
if show_norm:
- solver._solveBatchSize = 100
+ solver._solveBatchSize = 10
from plot_inf_set import plotInfSet1
- muInfVals, normEx, normApp, normErr = plotInfSet1(murange, murangeEff,
- approx, mu0, 250, 2.)
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet1(
+ murange, murangeEff, approx, mu0, 200, 2.)
+
diff --git a/examples/2d/pod/plot_inf_set.py b/examples/2d/pod/plot_inf_set.py
index 05a3ddc..1e3b68c 100644
--- a/examples/2d/pod/plot_inf_set.py
+++ b/examples/2d/pod/plot_inf_set.py
@@ -1,148 +1,223 @@
import warnings
import numpy as np
from matplotlib import pyplot as plt
-def plotInfSet1FromData(mus, Z, T, E, murange, approx, mu0, exp = 2.):
+def plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0, exp = 2.):
if hasattr(approx, "mus"):
mu2x = approx.mus(0) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]]
mu1 = np.real(np.power(mus, exp))
ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
+ Rmin, Rmax = np.min(R), np.max(R)
Emin, Emax = np.min(E), np.max(E)
+ betamin, betamax = np.min(beta), np.max(beta)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Z)
plt.semilogy(mu1, T, '--')
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [ZTmin, ZTmax], 'b:')
plt.plot(mu2x, [ZTmin] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [ZTmin, ZTmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [ZTmin, ZTmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|u(mu)|, |u_app(mu)|")
plt.grid()
plt.show()
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ plt.semilogy(mu1, R)
+ for l_ in approx.trainedModel.getPoles():
+ plt.plot([np.real(l_ ** exp)] * 2, [Rmin, Rmax], 'b:')
+ plt.plot(mu2x, [Rmax] * len(mu2x), 'kx')
+ plt.plot([murangeExp[0][0]] * 2, [Rmin, Rmax], 'm:')
+ plt.plot([murangeExp[1][0]] * 2, [Rmin, Rmax], 'm:')
+ plt.xlim(mu1[0], mu1[-1])
+ plt.title("|res(mu)|")
+ plt.grid()
+ plt.show()
+
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, E)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|err(mu)|")
plt.grid()
plt.show()
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ plt.semilogy(mu1, beta)
+ for l_ in approx.trainedModel.getPoles():
+ plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
+ plt.plot(mu2x, [betamax] * len(mu2x), 'kx')
+ plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
+ plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
+ plt.xlim(mu1[0], mu1[-1])
+ plt.title("beta(mu)")
+ plt.grid()
+ plt.show()
+
def plotInfSet1(murange, murangeEff, approx, mu0, nSamples = 20, exp = 2.):
mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
nSamples)
mus = np.power(mu1, 1. / exp)
+ F = approx.normRHS(mus)
Z = approx.normHF(mus)
T = approx.normApprox(mus)
+ R = approx.normRes(mus) / F
E = approx.normErr(mus) / Z
- plotInfSet1FromData(mus, Z, T, E, murange, approx, mu0, exp)
- return mus, Z, T, E
+ beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
+ plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0, exp)
+ return mus, Z, T, R, E, beta
-def plotInfSet2FromData(mus, Ze, Te, Ee, murange, approx, mu0,
+def plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange, approx, mu0,
exps = [2., 2.], clip = -1):
if hasattr(approx, "mus"):
mu2x, mu2y = approx.mus(0) ** exps[0], approx.mus(1) ** exps[1]
else:
mu2x, mu2y = mu0[0] ** exps[0], mu0[1] ** exps[1]
murangeExp = [[murange[0][0] ** exps[0], murange[0][1] ** exps[1]],
[murange[1][0] ** exps[0], murange[1][1] ** exps[1]]]
mu1s = np.unique([m[0] for m in mus])
mu2s = np.unique([m[1] for m in mus])
mu1 = np.power(mu1s, exps[0])
mu2 = np.power(mu2s, exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
- Z = np.log(Ze)
- T = np.log(Te)
- E = np.log(Ee)
+ Z = np.log10(Ze)
+ T = np.log10(Te)
+ R = np.log10(Re)
+ E = np.log10(Ee)
+ beta = np.log10(betae)
ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
+ Rmin, Rmax = np.min(R), np.max(R)
Emin, Emax = np.min(E), np.max(E)
+ betamin, betamax = np.min(beta), np.max(beta)
if clip > 0:
ZTmax -= clip * (ZTmax - ZTmin)
cmap = plt.cm.bone
else:
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = UserWarning)
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, Z, cmap = cmap,
levels = np.linspace(ZTmin, ZTmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Z, [ZTmin])
plt.plot(mu2x, mu2y, 'kx')
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.colorbar(p)
- plt.title("log|u(mu)|")
+ plt.title("log10|u(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, T, cmap = cmap,
levels = np.linspace(ZTmin, ZTmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, T, [ZTmin])
plt.plot(mu2x, mu2y, 'kx')
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
- plt.title("log|u_app(mu)|")
+ plt.title("log10|u_app(mu)|")
+ plt.colorbar(p)
+ plt.grid()
+ plt.show()
+
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ p = plt.contourf(Mu1, Mu2, R, levels = np.linspace(Rmin, Rmax, 50))
+ plt.plot(mu2x, mu2y, 'kx')
+ plt.plot([murangeExp[0][0]] * 2,
+ [murangeExp[0][1], murangeExp[1][1]], 'm:')
+ plt.plot([murangeExp[0][0], murangeExp[1][0]],
+ [murangeExp[1][1]] * 2, 'm:')
+ plt.plot([murangeExp[1][0]] * 2,
+ [murangeExp[1][1], murangeExp[0][1]], 'm:')
+ plt.plot([murangeExp[1][0], murangeExp[0][0]],
+ [murangeExp[0][1]] * 2, 'm:')
+ plt.title("log10|res(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, E, levels = np.linspace(Emin, Emax, 50))
plt.plot(mu2x, mu2y, 'kx')
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
- plt.title("log|err(mu)|")
+ plt.title("log10|err(mu)|")
+ plt.colorbar(p)
+ plt.grid()
+ plt.show()
+
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ p = plt.contourf(Mu1, Mu2, beta,
+ levels = np.linspace(betamin, betamax, 50))
+ plt.plot(mu2x, mu2y, 'kx')
+ plt.plot([murangeExp[0][0]] * 2,
+ [murangeExp[0][1], murangeExp[1][1]], 'm:')
+ plt.plot([murangeExp[0][0], murangeExp[1][0]],
+ [murangeExp[1][1]] * 2, 'm:')
+ plt.plot([murangeExp[1][0]] * 2,
+ [murangeExp[1][1], murangeExp[0][1]], 'm:')
+ plt.plot([murangeExp[1][0], murangeExp[0][0]],
+ [murangeExp[0][1]] * 2, 'm:')
+ plt.title("log10(beta(mu))")
plt.colorbar(p)
plt.grid()
plt.show()
def plotInfSet2(murange, murangeEff, approx, mu0, nSamples = 200,
exps = [2., 2.], clip = -1):
mu1 = np.linspace(murangeEff[0][0] ** exps[0], murangeEff[1][0] ** exps[0],
nSamples)
mu2 = np.linspace(murangeEff[0][1] ** exps[1], murangeEff[1][1] ** exps[1],
nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
mus = [(m1, m2) for m2 in mu2s for m1 in mu1s]
+ Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
Ze = approx.normHF(mus).reshape((nSamples, nSamples))
Te = approx.normApprox(mus).reshape((nSamples, nSamples))
+ Re = approx.normRes(mus).reshape((nSamples, nSamples)) / Fe
Ee = approx.normErr(mus).reshape((nSamples, nSamples)) / Ze
- plotInfSet2FromData(mus, Ze, Te, Ee, murange, approx, mu0, exps, clip)
- return mus, Ze, Te, Ee
+ betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus).reshape(
+ (nSamples, nSamples))
+ plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange,
+ approx, mu0, exps, clip)
+ return mus, Ze, Te, Re, Ee, betae
diff --git a/examples/2d/pod/plot_zero_set.py b/examples/2d/pod/plot_zero_set.py
index 0b0f622..69b32ed 100644
--- a/examples/2d/pod/plot_zero_set.py
+++ b/examples/2d/pod/plot_zero_set.py
@@ -1,80 +1,80 @@
import warnings
import numpy as np
from matplotlib import pyplot as plt
def plotZeroSet1(murange, murangeEff, approx, mu0, nSamples = 200, exp = 2.):
if hasattr(approx, "mus"):
mu2x = approx.mus(0) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]]
mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
nSamples)
mus = np.power(mu1, 1. / exp)
mu1 = np.real(mu1)
Z = approx.trainedModel.getQVal(mus)
Zabs = np.abs(Z)
Zmin, Zmax = np.min(Zabs), np.max(Zabs)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Zabs)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Zmin, Zmax], 'b--')
plt.plot(mu2x, [Zmax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [Zmin, Zmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Zmin, Zmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
def plotZeroSet2(murange, murangeEff, approx, mu0, nSamples = 200,
exps = [2., 2.], clip = -1):
if hasattr(approx, "mus"):
mu2x, mu2y = approx.mus(0) ** exps[0], approx.mus(1) ** exps[1]
else:
mu2x, mu2y = mu0[0] ** exps[0], mu0[1] ** exps[1]
murangeExp = [[murange[0][0] ** exps[0], murange[0][1] ** exps[1]],
[murange[1][0] ** exps[0], murange[1][1] ** exps[1]]]
mu1 = np.linspace(murangeEff[0][0] ** exps[0], murangeEff[1][0] ** exps[0],
nSamples)
mu2 = np.linspace(murangeEff[0][1] ** exps[1], murangeEff[1][1] ** exps[1],
nSamples)
mu1s = np.power(mu1, 1. / exps[0])
mu2s = np.power(mu2, 1. / exps[1])
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
mus = [(m1, m2) for m2 in mu2s for m1 in mu1s]
Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
- Zabs = np.log(np.abs(Z))
+ Zabs = np.log10(np.abs(Z))
Zabsmin, Zabsmax = np.min(Zabs), np.max(Zabs)
if clip > 0:
Zabsmin += clip * (Zabsmax - Zabsmin)
cmap = plt.cm.bone_r
else:
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = UserWarning)
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap,
levels = np.linspace(Zabsmin, Zabsmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Zabs, [Zabsmin])
plt.contour(Mu1, Mu2, np.real(Z), [0.], linestyles = 'dashed')
plt.contour(Mu1, Mu2, np.imag(Z), [0.], linewidths = 1,
linestyles = 'dotted')
plt.plot(mu2x, mu2y, 'kx')
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.colorbar(p)
- plt.title("log|Q(mu)|")
+ plt.title("log10|Q(mu)|")
plt.grid()
plt.show()
return mus, Z
diff --git a/examples/2d/pod/square_pod.py b/examples/2d/pod/square_pod.py
index ec45bca..c9768d0 100644
--- a/examples/2d/pod/square_pod.py
+++ b/examples/2d/pod/square_pod.py
@@ -1,150 +1,152 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = False
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
MN = 6
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [6.5 ** .5, .9 ** .5]
murange = [[6 ** .5, .5 ** .5], [8 ** .5, 1. ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
if ignore_forcing: solver.nbs = 1
rescaling = [lambda x: np.power(x, 2.), lambda x: x]
rescalingInv = [lambda x: np.power(x, .5), lambda x: x]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
- plotZeroSet2(murange, murangeEff, approx, mu0,
- 200, [2., 1.], clip = clip)
+ muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
+ 200, [2., 1.])
if show_norm:
+ solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- plotInfSet2(murange, murangeEff, approx, mu0,
- 25, [2., 1.], clip = clip)
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 25, [2., 1.], clip = clip)
+
diff --git a/examples/2d/pod/square_pod_hermite.py b/examples/2d/pod/square_pod_hermite.py
index 9454294..f512b92 100644
--- a/examples/2d/pod/square_pod_hermite.py
+++ b/examples/2d/pod/square_pod_hermite.py
@@ -1,105 +1,109 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
RandomSampler as RS,
ManualSampler as MS)
verb = 0
size = 1
show_sample = False
show_norm = True
MN = 4
R = (MN + 2) * (MN + 1) // 2
S0 = [3] * 2
S = [25]
assert R < np.prod(S)
samples = "centered"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [4.5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # large
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescaling = [lambda x: np.power(x, 2.), lambda x: x]
rescalingInv = [lambda x: np.power(x, .5), lambda x: x]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
method = RI
else:
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
else:
method = RBC
if samples == "distributed":
if sampling == "quadrature":
sampler0 = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
else: # if sampling == "random":
sampler0 = RS(murange, "SOBOL", scaling = rescaling,
scalingInv = rescalingInv)
S0 = np.prod(S0)
params['sampler'] = MS(murange, points = sampler0.generatePoints(S0))
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
- plotZeroSet2(murange, murangeEff, approx, mu0, 200, [2., 2.])
+ muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
+ 200, [2., 2.])
if show_norm:
+ solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- plotInfSet2(murange, murangeEff, approx, mu0, 25, [2., 2.])
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 25, [2., 2.])
+
diff --git a/examples/2d/pod/square_simplified_pod.py b/examples/2d/pod/square_simplified_pod.py
index 0af6d2f..41ac4d1 100644
--- a/examples/2d/pod/square_simplified_pod.py
+++ b/examples/2d/pod/square_simplified_pod.py
@@ -1,142 +1,146 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = False
show_norm = True
MN = 5
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # fat
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
elif size == 5: # tall
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]]
elif size == 6: # taller
mu0 = [11 ** .5, 2.25 ** .5]
mutar = [10.5 ** .5, 2.5 ** .5]
murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]]
elif size == 7: # low
mu0 = [7 ** .5, .75 ** .5]
mutar = [8 ** .5, 1 ** .5]
murange = [[6 ** .5, .25 ** .5], [8 ** .5, 1.25 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescaling = [lambda x: np.power(x, 2.)] * 2
rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
params['polybasis'] = "LEGENDRE"
params['polybasis'] = "MONOMIAL"
params['E'] = MN
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
- plotZeroSet2(murange, murangeEff, approx, mu0, 200, [2., 2.])
-
+ muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
+ 200, [2., 2.])
+
if show_norm:
+ solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- plotInfSet2(murange, murangeEff, approx, mu0, 25, [2., 2.])
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 25, [2., 2.])
+
diff --git a/examples/2d/pod/square_simplified_pod_hermite.py b/examples/2d/pod/square_simplified_pod_hermite.py
index 90e1593..d7bf757 100644
--- a/examples/2d/pod/square_simplified_pod_hermite.py
+++ b/examples/2d/pod/square_simplified_pod_hermite.py
@@ -1,105 +1,109 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
RandomSampler as RS,
ManualSampler as MS)
verb = 0
size = 1
show_sample = False
show_norm = True
MN = 4
R = (MN + 2) * (MN + 1) // 2
S0 = [3] * 2
S = [25]
assert R < np.prod(S)
samples = "centered"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "random"
if size == 1: # small
mu0 = [4 ** .5, 1.5 ** .5]
mutar = [5 ** .5, 1.75 ** .5]
murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]]
elif size == 2: # medium
mu0 = [4 ** .5, 1.75 ** .5]
mutar = [4.5 ** .5, 1.25 ** .5]
murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]]
elif size == 3: # large
mu0 = [6 ** .5, 4 ** .5]
mutar = [2 ** .5, 2.5 ** .5]
murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]]
elif size == 4: # crowded
mu0 = [10 ** .5, 2 ** .5]
mutar = [9 ** .5, 2.25 ** .5]
murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]]
aEff = 1.25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
(aEff*murange[0][1]**2. + bEff*murange[1][1]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
(aEff*murange[1][1]**2. + bEff*murange[0][1]**2.) ** .5]]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 20,
verbosity = verb)
rescaling = [lambda x: np.power(x, 2.)] * 2
rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
method = RI
else:
method = RP
else: #if algo == "RB":
params = {'R':R, 'S':S, 'POD':True}
if samples == "distributed":
method = RBD
else:
method = RBC
if samples == "distributed":
if sampling == "quadrature":
sampler0 = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
else: # if sampling == "random":
sampler0 = RS(murange, "SOBOL", scaling = rescaling,
scalingInv = rescalingInv)
S0 = np.prod(S0)
params['sampler'] = MS(murange, points = sampler0.generatePoints(S0))
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
if show_sample:
approx.plotApprox(mutar, name = 'u_app')
approx.plotHF(mutar, name = 'u_HF')
approx.plotErr(mutar, name = 'err')
approx.plotRes(mutar, name = 'res')
appErr, solNorm = approx.normErr(mutar), approx.normHF(mutar)
resNorm, RHSNorm = approx.normRes(mutar), approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational":
from plot_zero_set import plotZeroSet2
- plotZeroSet2(murange, murangeEff, approx, mu0, 200, [2., 2.])
+ muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
+ 200, [2., 2.])
if show_norm:
+ solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- plotInfSet2(murange, murangeEff, approx, mu0, 25, [2., 2.])
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 50, [2., 2.])
+
diff --git a/examples/2d/pod/cookie_single_pod.py b/examples/2d/pod/synthetic_pod.py
similarity index 76%
copy from examples/2d/pod/cookie_single_pod.py
copy to examples/2d/pod/synthetic_pod.py
index 752762f..ce5ef62 100644
--- a/examples/2d/pod/cookie_single_pod.py
+++ b/examples/2d/pod/synthetic_pod.py
@@ -1,138 +1,157 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
- CookieEngineSingle as CES
+ SyntheticBivariateEngine as SBE
from rrompy.reduction_methods.centered import RationalPade as RP
from rrompy.reduction_methods.distributed import RationalInterpolant as RI
from rrompy.reduction_methods.centered import RBCentered as RBC
from rrompy.reduction_methods.distributed import RBDistributed as RBD
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
clip = -1
#clip = .4
#clip = .6
-MN = 15
+Delta = 0
+MN = 10
R = (MN + 2) * (MN + 1) // 2
S = [int(np.ceil(R ** .5))] * 2
+PODTol = 1e-10
samples = "centered"
samples = "centered_fake"
samples = "distributed"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
-sampling = "quadrature_total"
-sampling = "random"
+#sampling = "quadrature_total"
+#sampling = "random"
-if size == 1: # below
- mu0 = [20 ** .5, 1. ** .5]
- mutar = [20.5 ** .5, 1.05 ** .5]
- murange = [[18.5 ** .5, .85 ** .5], [21.5 ** .5, 1.15 ** .5]]
+if samples == "distributed":
+ radial = 0
+# radial = "gaussian"
+# radial = "thinplate"
+# radial = "multiquadric"
+ rW0 = 10.
+ radialWeight = [r * rW0 for r in [5., 5.]]
+
+if size == 1: # small
+ mu0 = [10. ** .5, 15. ** .5]
+ mutar = [12. ** .5, 14. ** .5]
+ murange = [[5. ** .5, 10. ** .5], [15 ** .5, 20 ** .5]]
+if size == 2: # large
+ mu0 = [15. ** .5, 17.5 ** .5]
+ mutar = [18. ** .5, 22. ** .5]
+ murange = [[5. ** .5, 10. ** .5], [25 ** .5, 25 ** .5]]
+if size == 3: # medium
+ mu0 = [17.5 ** .5, 15 ** .5]
+ mutar = [20. ** .5, 18. ** .5]
+ murange = [[10. ** .5, 10. ** .5], [25 ** .5, 20 ** .5]]
+
+assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
kappa = 20. ** .5
theta = - np.pi / 6.
-n = 40
-Rad = 1.
+n = 20
L = np.pi
-nX = 2
-nY = 1
-solver = CES(kappa = kappa, theta = theta, n = n, R = Rad, L = L, nX = nX,
- nY = nY, mu0 = mu0, verbosity = verb)
+solver = SBE(kappa = kappa, theta = theta, n = n, L = L,
+ mu0 = mu0, verbosity = verb)
rescaling = [lambda x: np.power(x, 2.)] * 2
rescalingInv = [lambda x: np.power(x, .5)] * 2
if algo == "rational":
- params = {'N':MN, 'M':MN, 'S':S, 'POD':True}
+ params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
if samples == "distributed":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
params['E'] = MN
+ params['radialBasis'] = radial
+ params['radialBasisWeights'] = radialWeight
method = RI
elif samples == "centered_fake":
params['polybasis'] = "MONOMIAL"
params['S'] = R
method = RI
else:
params['S'] = R
method = RP
else: #if algo == "RB":
- params = {'R':R, 'S':S, 'POD':True}
+ params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S,
+ 'POD':True, 'PODTolerance':PODTol}
if samples == "distributed":
method = RBD
elif samples == "centered_fake":
params['S'] = R
method = RBD
else:
params['S'] = R
method = RBC
if samples == "distributed":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
# scalingInv = rescalingInv)
# params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
# scalingInv = rescalingInv)
params['S'] = [max(j, MN + 1) for j in params['S']]
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
scalingInv = rescalingInv)
params['S'] = R
elif samples == "centered_fake":
params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
scalingInv = rescalingInv)
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', homogeneized = False,
what = "REAL")
approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL")
approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL")
# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
if algo == "rational" and approx.N > 0:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 2.], clip = clip)
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
- muInfVals, normEx, normApp, normErr = plotInfSet2(murange, murangeEff,
- approx, mu0, 40,
- [2., 2.], clip = clip)
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 50, [2., 2.], clip = clip)
diff --git a/rrompy/hfengines/base/matrix_engine_base.py b/rrompy/hfengines/base/matrix_engine_base.py
index 5908ff3..ae63291 100644
--- a/rrompy/hfengines/base/matrix_engine_base.py
+++ b/rrompy/hfengines/base/matrix_engine_base.py
@@ -1,460 +1,498 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
import scipy.sparse as scsp
from matplotlib import pyplot as plt
from copy import deepcopy as copy, copy as softcopy
from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, strLst, Tuple, List,
DictAny, paramVal, paramList,
sampList)
from rrompy.utilities.base import (purgeList, getNewFilename, verbosityDepth,
multibinom)
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList, emptySampleList
from rrompy.solver import setupSolver
from rrompy.solver.scipy import tensorizeLS, detensorizeLS
__all__ = ['MatrixEngineBase']
class MatrixEngineBase:
"""
Generic solver for parametric matrix problems.
Attributes:
verbosity: Verbosity level.
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.
"""
_solveBatchSize = 1
def __init__(self, verbosity : int = 10, timestamp : bool = True):
self.verbosity = verbosity
self.timestamp = timestamp
self.nAs, self.nbs = 1, 1
self.setSolver("SPSOLVE", {"use_umfpack" : False})
self.npar = 0
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 __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def __deepcopy__(self, memo):
return softcopy(self)
@property
def npar(self):
"""Value of npar."""
return self._npar
@npar.setter
def npar(self, npar):
nparOld = self._npar if hasattr(self, "_npar") else -1
if npar != nparOld:
self.rescalingExp = [1.] * npar
self._npar = npar
@property
def nAs(self):
"""Value of nAs."""
return self._nAs
@nAs.setter
def nAs(self, nAs):
self._nAs = nAs
self.resetAs()
@property
def nbs(self):
"""Value of nbs."""
return self._nbs
@nbs.setter
def nbs(self, nbs):
self._nbs = nbs
self.resetbs()
@property
def nbsH(self) -> int:
return max(self.nbs, self.nAs)
def spacedim(self):
return self.As[0].shape[1]
def checkParameter(self, mu:paramList):
return checkParameter(mu, self.npar)
def checkParameterList(self, mu:paramList):
return checkParameterList(mu, self.npar)
def buildEnergyNormForm(self): # eye
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = np.eye(self.spacedim())
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D:
"""Scalar product."""
if not hasattr(self, "energyNormMatrix"):
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling energy matrix.",
timestamp = self.timestamp)
self.buildEnergyNormForm()
if self.verbosity >= 20:
verbosityDepth("DEL", "Done assembling energy matrix.",
timestamp = self.timestamp)
if not isinstance(u, (np.ndarray,)): u = u.data
if not isinstance(v, (np.ndarray,)): v = v.data
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 checkAInBounds(self, derI : int = 0):
"""Check if derivative index is oob for operator of linear system."""
if derI < 0 or derI >= self.nAs:
d = self.spacedim()
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
def checkbInBounds(self, derI : int = 0, homogeneized : bool = False):
"""Check if derivative index is oob for RHS of linear system."""
nbs = self.nbsH if homogeneized else self.nbs
if derI < 0 or derI >= nbs:
return np.zeros(self.spacedim(), dtype = np.complex)
def resetAs(self):
"""Reset (derivatives of) operator of linear system."""
self.setAs([None] * self.nAs)
if hasattr(self, "_nbs"): self.resetbsH()
def resetbs(self):
"""Reset (derivatives of) RHS of linear system."""
self.setbs([None] * self.nbs)
if hasattr(self, "_nAs"): self.resetbsH()
def resetbsH(self):
"""Reset (derivatives of) homogeneized RHS of linear system."""
self.setbsH([None] * self.nbsH)
def setAs(self, As:List[Np2D]):
"""Assign terms of operator of linear system."""
if len(As) != self.nAs:
raise RROMPyException(("Expected number {} of terms of As not "
"matching given list length {}.").format(self.nAs,
len(As)))
self.As = [copy(A) for A in As]
def setbs(self, bs:List[Np1D]):
"""Assign terms of RHS of linear system."""
if len(bs) != self.nbs:
raise RROMPyException(("Expected number {} of terms of bs not "
"matching given list length {}.").format(self.nbs,
len(bs)))
self.bs = [copy(b) for b in bs]
def setbsH(self, bsH:List[Np1D]):
"""Assign terms of homogeneized RHS of linear system."""
if len(bsH) != self.nbsH:
raise RROMPyException(("Expected number {} of terms of bsH not "
"matching given list length {}.").format(self.nbsH,
len(bsH)))
self.bsH = [copy(bH) for bH in bsH]
def _assembleA(self, mu : paramVal = [], der : List[int] = 0,
derI : int = None, muBase : paramVal = None) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if muBase is None: muBase = [0.] * self.npar
muBase = self.checkParameter(muBase)
if self.npar > 0: mu, muBase = mu[0], muBase[0]
if not hasattr(der, "__len__"): der = [der] * self.npar
if derI is None: derI = hashD(der)
Anull = self.checkAInBounds(derI)
if Anull is not None: return Anull
rExp = self.rescalingExp
A = copy(self.As[derI])
for j in range(derI + 1, self.nAs):
derIdx = hashI(j, self.npar)
diffIdx = [x - y for (x, y) in zip(derIdx, der)]
if np.all([x >= 0 for x in diffIdx]):
A = A + (np.prod((mu ** rExp - muBase ** rExp) ** diffIdx)
* multibinom(derIdx, diffIdx) * self.As[j])
return A
@abstractmethod
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""
Assemble terms of operator of linear system and return it (or its
derivative) at a given parameter.
"""
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
for j in range(derI, self.nAs):
if self.As[j] is None: self.As[j] = self.checkAInBounds(-1)
return self._assembleA(mu, der, derI)
def affineLinearSystemA(self, mu : paramVal = []) -> List[Np2D]:
"""
Assemble affine blocks of operator of linear system (just linear
blocks).
"""
As = [None] * self.nAs
for j in range(self.nAs):
As[j] = self.A(mu, hashI(j, self.npar))
return As
def affineWeightsA(self, mu : paramVal = []) -> List[str]:
"""
Assemble affine blocks of operator of linear system (just affine
weights). Stored as strings for the sake of pickling.
"""
mu = self.checkParameter(mu)
lambdasA = ["1."]
mu0Eff = mu ** self.rescalingExp
for j in range(1, self.nAs):
lambdasA += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format(
','.join([str(x) for x in mu0Eff[0]]),
self.rescalingExp, hashI(j, self.npar))]
return lambdasA
def affineBlocksA(self, mu : paramVal = [])\
-> Tuple[List[Np2D], List[str]]:
"""Assemble affine blocks of operator of linear system."""
return self.affineLinearSystemA(mu), self.affineWeightsA(mu)
def _assembleb(self, mu : paramVal = [], der : List[int] = 0,
derI : int = None, homogeneized : bool = False,
muBase : paramVal = None) -> ScOp:
"""Assemble (derivative of) (homogeneized) RHS of linear system."""
mu = self.checkParameter(mu)
if muBase is None: muBase = [0.] * self.npar
muBase = self.checkParameter(muBase)
if self.npar > 0: mu, muBase = mu[0], muBase[0]
if not hasattr(der, "__len__"): der = [der] * self.npar
if derI is None: derI = hashD(der)
bnull = self.checkbInBounds(derI, homogeneized)
if bnull is not None: return bnull
bs = self.bsH if homogeneized else self.bs
rExp = self.rescalingExp
b = copy(bs[derI])
for j in range(derI + 1, len(bs)):
derIdx = hashI(j, self.npar)
diffIdx = [x - y for (x, y) in zip(derIdx, der)]
if np.all([x >= 0 for x in diffIdx]):
b = b + (np.prod((mu ** rExp - muBase ** rExp) ** diffIdx)
* multibinom(derIdx, diffIdx) * bs[j])
return b
@abstractmethod
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""
Assemble terms of (homogeneized) RHS of linear system and return it (or
its derivative) at a given parameter.
"""
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
if homogeneized:
for j in range(derI, self.nbsH):
if self.bsH[j] is None: self.bsH[j] = self.checkbInBounds(-1)
else:
for j in range(derI, self.nbs):
if self.bs[j] is None: self.bs[j] = self.checkbInBounds(-1)
return self._assembleb(mu, der, derI, homogeneized)
def affineLinearSystemb(self, mu : paramVal = [],
homogeneized : bool = False) -> List[Np1D]:
"""
Assemble affine blocks of RHS of linear system (just linear blocks).
"""
nbs = self.nbsH if homogeneized else self.nbs
bs = [None] * nbs
for j in range(nbs):
bs[j] = self.b(mu, hashI(j, self.npar), homogeneized)
return bs
def affineWeightsb(self, mu : paramVal = [],
homogeneized : bool = False) -> List[str]:
"""
Assemble affine blocks of RHS of linear system (just affine weights).
Stored as strings for the sake of pickling.
"""
mu = self.checkParameter(mu)
nbs = self.nbsH if homogeneized else self.nbs
lambdasb = ["1."]
mu0Eff = mu ** self.rescalingExp
for j in range(1, nbs):
lambdasb += ["np.prod((mu ** ({1}) - [{0}]) ** ({2}))".format(
','.join([str(x) for x in mu0Eff[0]]),
self.rescalingExp, hashI(j, self.npar))]
return lambdasb
def affineBlocksb(self, mu : paramVal = [], homogeneized : bool = False)\
-> Tuple[List[Np1D], List[str]]:
"""Assemble affine blocks of RHS of linear system."""
return (self.affineLinearSystemb(mu, homogeneized),
self.affineWeightsb(mu, homogeneized))
def setSolver(self, solverType:str, solverArgs : DictAny = {}):
"""Choose solver type and parameters."""
self._solver, self._solverArgs = setupSolver(solverType, solverArgs)
def solve(self, mu : paramList = [], RHS : sampList = None,
homogeneized : bool = False) -> sampList:
"""
Find solution of linear system.
Args:
mu: parameter value.
RHS: RHS of linear system. If None, defaults to that of parametric
system. Defaults to None.
"""
mu, _ = self.checkParameterList(mu)
if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype)
sol = emptySampleList()
if len(mu) > 0:
if RHS is None:
RHS = [self.b(m, homogeneized = homogeneized) for m in mu]
RHS = sampleList(RHS)
mult = 0 if len(RHS) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size")
u = self._solver(self.A(mu[0]), RHS[0], self._solverArgs)
sol.reset((len(u), len(mu)), dtype = u.dtype)
sol[0] = u
for j in range(1, len(mu), self._solveBatchSize):
if self._solveBatchSize != 1:
iRange = list(range(j, min(j + self._solveBatchSize,
len(mu))))
As = [self.A(mu[i]) for i in iRange]
bs = [RHS[mult * i] for i in iRange]
A, b = tensorizeLS(As, bs)
else:
A, b = self.A(mu[j]), RHS[mult * j]
solStack = self._solver(A, b, self._solverArgs)
if self._solveBatchSize != 1:
sol[iRange] = detensorizeLS(solStack, len(iRange))
else:
sol[j] = solStack
return sol
def residual(self, u:sampList, mu : paramList = [],
homogeneized : bool = False) -> sampList:
"""
Find residual of linear system for given approximate solution.
Args:
u: numpy complex array with function dofs. If None, set to 0.
mu: parameter value.
"""
mu, _ = self.checkParameterList(mu)
if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype)
if u is not None:
u = sampleList(u)
mult = 0 if len(u) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size")
res = emptySampleList()
for j in range(len(mu)):
b = self.b(mu[j], homogeneized = homogeneized)
if u is None:
r = b
else:
r = b - self.A(mu[j]).dot(u[mult * j])
if j == 0:
res.reset((len(r), len(mu)), dtype = r.dtype)
res[j] = r
return res
+ def _rayleighQuotient(self, A:Np2D, v0:Np1D, M:Np2D, sigma : float = 0.,
+ nIterP : int = 10, nIterR : int = 10) -> float:
+ v0 /= v0.T.conj().dot(M.dot(v0)) ** .5
+ for j in range(nIterP):
+ v0 = self._solver(A - sigma * M, M.dot(v0), self._solverArgs)
+ v0 /= v0.T.conj().dot(M.dot(v0)) ** .5
+ l0 = v0.T.conj().dot(A.dot(v0))
+ for j in range(nIterR):
+ v0 = self._solver(A - l0 * M, M.dot(v0), self._solverArgs)
+ v0 /= v0.T.conj().dot(M.dot(v0)) ** .5
+ l0 = v0.T.conj().dot(A.dot(v0))
+ if np.isnan(l0): l0 = 0.
+ return np.abs(l0)
+
+ def stabilityFactor(self, u:sampList, mu : paramList = [],
+ nIterP : int = 10, nIterR : int = 10) -> sampList:
+ """
+ Find stability factor of matrix of linear system using iterative
+ inverse power iteration- and Rayleigh quotient-based procedure.
+
+ Args:
+ u: numpy complex arrays with function dofs.
+ mu: parameter values.
+ nIterP: number of iterations of power method.
+ nIterR: number of iterations of Rayleigh quotient method.
+ """
+ mu, _ = self.checkParameterList(mu)
+ if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype)
+ u = sampleList(u)
+ mult = 0 if len(u) == 1 else 1
+ RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size")
+ stabFact = np.empty(len(mu), dtype = float)
+ for j in range(len(mu)):
+ stabFact[j] = self._rayleighQuotient(self.A(mu[j]), u[mult * j],
+ self.energyNormMatrix,
+ 0., nIterP, nIterR)
+ return stabFact
+
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
idxs = np.arange(self.spacedim())
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.abs(u).flatten())
plt.title("|{0}|".format(name))
if 'PHASE' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.angle(u).flatten())
plt.title("phase({0})".format(name))
if 'REAL' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.real(u).flatten())
plt.title("Re({0})".format(name))
if 'IMAG' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.imag(u).flatten())
plt.title("Im({0})".format(name))
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()
diff --git a/rrompy/hfengines/linear_problem/bidimensional/__init__.py b/rrompy/hfengines/linear_problem/bidimensional/__init__.py
index abf4add..54ca3ff 100644
--- a/rrompy/hfengines/linear_problem/bidimensional/__init__.py
+++ b/rrompy/hfengines/linear_problem/bidimensional/__init__.py
@@ -1,36 +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 .cookie_engine_single import CookieEngineSingle
from .membrane_fracture_engine import MembraneFractureEngine
from .laplace_disk_gaussian_2 import LaplaceDiskGaussian2
from .helmholtz_square_simplified_domain_problem_engine import \
HelmholtzSquareSimplifiedDomainProblemEngine
from .helmholtz_square_domain_problem_engine import \
HelmholtzSquareDomainProblemEngine
+from .synthetic_bivariate_engine import SyntheticBivariateEngine
__all__ = [
'CookieEngineSingle',
'MembraneFractureEngine',
'LaplaceDiskGaussian2',
'HelmholtzSquareSimplifiedDomainProblemEngine',
- 'HelmholtzSquareDomainProblemEngine'
+ 'HelmholtzSquareDomainProblemEngine',
+ 'SyntheticBivariateEngine'
]
diff --git a/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py b/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
new file mode 100644
index 0000000..21f2b87
--- /dev/null
+++ b/rrompy/hfengines/linear_problem/bidimensional/synthetic_bivariate_engine.py
@@ -0,0 +1,153 @@
+# 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
+import ufl
+from rrompy.utilities.base.types import ScOp, List, paramVal
+from rrompy.solver.fenics import fenONE, fenZERO
+from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
+ HelmholtzProblemEngine)
+from rrompy.utilities.base import verbosityDepth
+from rrompy.utilities.poly_fitting.polynomial import (
+ hashDerivativeToIdx as hashD)
+
+__all__ = ['SyntheticBivariateEngine']
+
+class SyntheticBivariateEngine(HelmholtzProblemEngine):
+
+ def __init__(self, kappa:float, theta:float, n:int, L : int = 2.,
+ mu0 : paramVal = [12. ** .5, 15. ** .5],
+ degree_threshold : int = np.inf, verbosity : int = 10,
+ timestamp : bool = True):
+ super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
+ verbosity = verbosity, timestamp = timestamp)
+ self.nAs = 3
+ self.npar = 2
+ self.rescalingExp = [2., 2.]
+
+ mesh = fen.RectangleMesh(fen.Point(0., 0.), fen.Point(L, L), n, n)
+ self.V = fen.FunctionSpace(mesh, "P", 1)
+
+ x, y = fen.SpatialCoordinate(mesh)[:]
+ self._above = ufl.conditional(ufl.ge(y, .5 * L), fenONE, fenZERO)
+ self._below = fenONE - self._above
+ c, s = np.cos(theta), np.sin(theta)
+ self.forcingTerm = [fen.cos(kappa * (c * x + s * y)),
+ fen.sin(kappa * (c * x + s * y))]
+
+ def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
+ """Assemble (derivative of) operator of linear system."""
+ mu = self.checkParameter(mu)
+ if not hasattr(der, "__len__"): der = [der] * self.npar
+ derI = hashD(der)
+ self.autoSetDS()
+ if derI <= 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 derI <= 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 * self._above * fen.dot(self.u, self.v) * fen.dx
+ a1Im = - n2Im * self._above * 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 derI <= 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 * self._below * fen.dot(self.u, self.v) * fen.dx
+ a2Im = - n2Im * self._below * 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)
+ return self._assembleA(mu, der, derI)