diff --git a/examples/2d/base/fracture.py b/examples/2d/base/fracture.py
index af49ae8..b98b5e5 100644
--- a/examples/2d/base/fracture.py
+++ b/examples/2d/base/fracture.py
@@ -1,46 +1,46 @@
import numpy as np
import ufl
import fenics as fen
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.solver.fenics import affine_warping
verb = 100
mu0 = [45. ** .5, .7]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta,
n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu0)[0]
#solver.plotmesh(figsize = (7.5, 4.5))
#solver.plot(u0, what = 'REAL', figsize = (8, 5))
print(solver.norm(uh))
#solver.plot(uh, what = 'REAL', figsize = (8, 5))
-#solver.plot(solver.residual(uh, mu0)[0], name = 'res',
+#solver.plot(solver.residual(mu0, uh)[0], name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], filename = 'out', folder = True,
# forceNewFile = False)
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * mu0[1]]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * mu0[1]]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
#solver.plotmesh([warp, warpI], figsize = (7.5, 4.5))
#solver.plot(u0, [warp, warpI], what = 'REAL', figsize = (8, 5))
solver.plot(uh, [warp, warpI], what = 'REAL', figsize = (8, 5))
-#solver.plot(solver.residual(uh, mu0)[0], [warp, warpI], name = 'res',
+#solver.plot(solver.residual(mu0, uh)[0], [warp, warpI], name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], [warp, warpI],
# filename = 'outW', folder = True,
# forceNewFile = False)
diff --git a/examples/2d/base/fracture_nodomain.py b/examples/2d/base/fracture_nodomain.py
index 3163f30..400c730 100644
--- a/examples/2d/base/fracture_nodomain.py
+++ b/examples/2d/base/fracture_nodomain.py
@@ -1,47 +1,47 @@
import numpy as np
import ufl
import fenics as fen
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.solver.fenics import affine_warping
verb = 100
mu0Aug = [45. ** .5, .6]
mu0Aug = [45. ** .5, .1]
mu0 = mu0Aug[0]
H = 1.
L = .75
delta = .05
n = 50
solver = MFEND(mu0 = mu0Aug, H = H, L = L, delta = delta,
n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu0)[0]
solver.plotmesh(figsize = (7.5, 4.5))
solver.plot(u0, what = 'REAL', figsize = (8, 5))
print(solver.norm(uh))
solver.plot(uh, what = 'REAL', figsize = (8, 5))
-solver.plot(solver.residual(uh, mu0)[0], name = 'res',
+solver.plot(solver.residual(mu0, uh)[0], name = 'res',
what = 'REAL', figsize = (8, 5))
solver.outParaviewTimeDomain(uh, mu0, filename = 'outND', folder = True)
##
L = mu0Aug[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)
solver.plotmesh([warp, warpI], figsize = (7.5, 4.5))
solver.plot(u0, [warp, warpI], what = 'REAL', figsize = (8, 5))
solver.plot(uh, [warp, warpI], what = 'REAL', figsize = (8, 5))
-solver.plot(solver.residual(uh, mu0)[0], [warp, warpI], name = 'res',
+solver.plot(solver.residual(mu0, uh)[0], [warp, warpI], name = 'res',
what = 'REAL', figsize = (8, 5))
solver.outParaviewTimeDomain(uh, mu0, [warp, warpI],
filename = 'outNDW', folder = True)
diff --git a/examples/2d/base/square.py b/examples/2d/base/square.py
index ae2428e..856cf1a 100644
--- a/examples/2d/base/square.py
+++ b/examples/2d/base/square.py
@@ -1,13 +1,13 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
verb = 100
mu0 = [4 ** .5, 1.5 ** .5]
solver = HSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 50,
verbosity = verb)
uh = solver.solve(mu0)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
-solver.plot(solver.residual(uh, mu0)[0], 'res')
+solver.plot(solver.residual(mu0, uh)[0], 'res')
diff --git a/examples/2d/base/square_simplified.py b/examples/2d/base/square_simplified.py
index fa0e984..6188ee7 100644
--- a/examples/2d/base/square_simplified.py
+++ b/examples/2d/base/square_simplified.py
@@ -1,13 +1,13 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareSimplifiedDomainProblemEngine as HSSDPE
verb = 0
mu0 = [4 ** .5, 1.5 ** .5]
solver = HSSDPE(kappa = 2.5, theta = np.pi / 3, mu0 = mu0, n = 50,
verbosity = verb)
uh = solver.solve(mu0)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
-solver.plot(solver.residual(uh, mu0)[0], 'res')
+solver.plot(solver.residual(mu0, uh)[0], 'res')
diff --git a/examples/2d/greedy/fracture_pod.py b/examples/2d/greedy/fracture_pod.py
new file mode 100644
index 0000000..ea4eb71
--- /dev/null
+++ b/examples/2d/greedy/fracture_pod.py
@@ -0,0 +1,223 @@
+import numpy as np
+from rrompy.hfengines.linear_problem.bidimensional import \
+ MembraneFractureEngine as MFE
+from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RIG
+from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RBG
+from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
+ QuadratureSamplerTotal as QST)
+
+verb = 15
+size = 2
+show_sample = True
+show_norm = True
+
+MN = 1
+S = (MN + 2) * (MN + 1) // 2
+greedyTol = 1e-2
+collinearityTol = 0.
+nTestPoints = 400
+
+algo = "rational"
+algo = "RB"
+
+if algo == "rational":
+ radial = ""
+ #radial = "_gaussian"
+ #radial = "_thinplate"
+ #radial = "_multiquadric"
+ rW0 = 5.
+ radialWeight = [rW0] * 2
+ polybasis = "CHEBYSHEV"
+ polybasis = "LEGENDRE"
+ #polybasis = "MONOMIAL"
+ errorEstimatorKind = 'AFFINE'
+ errorEstimatorKind = 'DISCREPANCY'
+# errorEstimatorKind = 'INTERPOLATORY'
+# errorEstimatorKind = 'LOCALL2'
+# errorEstimatorKind = 'DIAGONAL'
+
+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 == 11: #L below
+ mu0 = [110 ** .5, .4]
+ mutar = [115 ** .5, .4]
+ murange = [[90 ** .5, .3], [130 ** .5, .5]]
+elif size == 12: #L top
+ mu0 = [110 ** .5, .6]
+ mutar = [115 ** .5, .6]
+ murange = [[90 ** .5, .5], [130 ** .5, .7]]
+elif size == 13: #L interesting
+ mu0 = [110 ** .5, .5]
+ mutar = [115 ** .5, .5]
+ murange = [[90 ** .5, .3], [130 ** .5, .7]]
+elif size == 14: #L belowbelow
+ mu0 = [110 ** .5, .2]
+ mutar = [115 ** .5, .2]
+ murange = [[90 ** .5, .1], [130 ** .5, .3]]
+elif size == 15: #L toptop
+ mu0 = [110 ** .5, .8]
+ mutar = [115 ** .5, .8]
+ murange = [[90 ** .5, .7], [130 ** .5, .9]]
+elif size == 16: #L interestinginteresting
+ mu0 = [110 ** .5, .5]
+ mutar = [115 ** .5, .6]
+ murange = [[90 ** .5, .1], [130 ** .5, .9]]
+elif size == 17: #L interestingtop
+ mu0 = [110 ** .5, .7]
+ mutar = [115 ** .5, .6]
+ murange = [[90 ** .5, .5], [130 ** .5, .9]]
+elif size == 18: #L interestingbelow
+ mu0 = [110 ** .5, .3]
+ mutar = [115 ** .5, .4]
+ murange = [[90 ** .5, .1], [130 ** .5, .5]]
+elif size == 100: # tiny
+ mu0 = [32.5 ** .5, .5]
+ mutar = [34 ** .5, .5]
+ murange = [[30 ** .5, .3], [35 ** .5, .7]]
+
+aEff = 1#.05
+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 = 20
+
+solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
+
+rescalingExp = [2., 1.]
+if algo == "rational":
+ params = {'S':S, 'POD':True, 'greedyTol':greedyTol,
+ 'sampler':QS(murange, 'UNIFORM', scalingExp = rescalingExp),
+ 'nTestPoints':nTestPoints, 'polybasis':polybasis + radial,
+ 'radialDirectionalWeights':radialWeight,
+ 'errorEstimatorKind':errorEstimatorKind,
+ 'trainSetGenerator':QST(murange, 'CHEBYSHEV',
+ scalingExp = rescalingExp)}
+ method = RIG
+else: #if algo == "RB":
+ params = {'S':S, 'POD':True, 'greedyTol':greedyTol,
+ 'sampler':QS(murange, 'UNIFORM', scalingExp = rescalingExp),
+ 'nTestPoints':nTestPoints,
+ 'trainSetGenerator':QST(murange, 'CHEBYSHEV',
+ scalingExp = rescalingExp)}
+ method = RBG
+
+approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
+approx.samplingEngine.allowRepeatedSamples = False
+
+approx.greedy(True)
+
+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', what = "REAL")
+ approx.plotHF(mutar, [warp, warpI], name = 'u_HF', what = "REAL")
+ approx.plotErr(mutar, [warp, warpI], name = 'err', what = "REAL")
+# approx.plotRes(mutar, [warp, warpI], name = 'res', 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)))
+
+verb = approx.trainedModel.verbosity
+approx.trainedModel.verbosity = 5
+try:
+ from plot_zero_set import plotZeroSet2
+ muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
+ 200, [2., 1.])
+except: pass
+
+if show_norm:
+ solver._solveBatchSize = 100
+ from plot_inf_set import plotInfSet2
+ muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
+ murange, murangeEff, approx, mu0, 50,
+ [2., 1.], relative = True,
+ nobeta = True)
+
+ import warnings
+ from matplotlib import pyplot as plt
+ mu2x, mu2y = approx.mus(0) ** 2., approx.mus(1) ** 1.
+ murangeExp = [[murange[0][0] ** 2., murange[0][1]],
+ [murange[1][0] ** 2., murange[1][1]]]
+ mu1s = np.unique([m[0] for m in muInfVals])
+ mu2s = np.unique([m[1] for m in muInfVals])
+ mu1 = np.power(mu1s, 2.)
+ mu2 = np.power(mu2s, 1.)
+ Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
+ Reste = (approx.errorEstimator(muInfVals)
+ * approx.trainedModel.getQVal(muInfVals)).reshape(normEx.shape)
+ Rest = np.log10(Reste)
+ Restmin, Restmax = np.min(Rest), np.max(Rest)
+ cmap = plt.cm.jet
+ warnings.simplefilter("ignore", category = (UserWarning,
+ np.ComplexWarning))
+ plt.figure(figsize = (15, 7))
+ plt.jet()
+ p = plt.contourf(Mu1, Mu2, Rest,
+ levels = np.linspace(Restmin, Restmax, 50))
+ plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
+ for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
+ plt.annotate("{}".format(j), xy)
+ 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_est(mu)|")
+ plt.colorbar(p)
+ plt.grid()
+ plt.show()
+approx.trainedModel.verbosity = verb
+
+try:
+ print(np.sort(approx.getPoles([None, .5]) ** 2.))
+except: pass
diff --git a/examples/2d/greedy/plot_inf_set.py b/examples/2d/greedy/plot_inf_set.py
new file mode 120000
index 0000000..ad92d4d
--- /dev/null
+++ b/examples/2d/greedy/plot_inf_set.py
@@ -0,0 +1 @@
+/home/pradovera/Repos/RROMPy/examples/2d/pod/plot_inf_set.py
\ No newline at end of file
diff --git a/examples/2d/greedy/plot_zero_set.py b/examples/2d/greedy/plot_zero_set.py
new file mode 120000
index 0000000..6b5351e
--- /dev/null
+++ b/examples/2d/greedy/plot_zero_set.py
@@ -0,0 +1 @@
+/home/pradovera/Repos/RROMPy/examples/2d/pod/plot_zero_set.py
\ No newline at end of file
diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py
index 35463e0..3f102c2 100644
--- a/examples/2d/pod/fracture_pod.py
+++ b/examples/2d/pod/fracture_pod.py
@@ -1,261 +1,267 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
-from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
+#from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 10
size = 2
show_sample = True
show_norm = True
Delta = 0
MN = 5
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
-PODTol = 1e-6
+PODTol = 1e-10
SPivot = MN + 1
SMarginal = 4
MMarginal = SMarginal - 1
matchingWeight = 1.
cutOffTolerance = 1. * np.inf
cutOffType = "MAGNITUDE"
samples = "centered"
samples = "standard"
#samples = "pivoted"
algo = "rational"
-#algo = "RB"
+algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
polydegreetype = "TOTAL"
polydegreetype = "FULL"
if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
if samples == "pivoted":
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
assert Delta <= 0
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 == 11: #L below
mu0 = [110 ** .5, .4]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .3], [130 ** .5, .5]]
elif size == 12: #L top
mu0 = [110 ** .5, .6]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .7]]
elif size == 13: #L interesting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .5]
murange = [[90 ** .5, .3], [130 ** .5, .7]]
elif size == 14: #L belowbelow
mu0 = [110 ** .5, .2]
mutar = [115 ** .5, .2]
murange = [[90 ** .5, .1], [130 ** .5, .3]]
elif size == 15: #L toptop
mu0 = [110 ** .5, .8]
mutar = [115 ** .5, .8]
murange = [[90 ** .5, .7], [130 ** .5, .9]]
elif size == 16: #L interestinginteresting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .1], [130 ** .5, .9]]
elif size == 17: #L interestingtop
mu0 = [110 ** .5, .7]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .9]]
elif size == 18: #L interestingbelow
mu0 = [110 ** .5, .3]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .1], [130 ** .5, .5]]
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 = 20
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True,
'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight * 2
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "pivoted":
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RIP
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':R,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
method = RB
elif samples == "pivoted":
- params['S'] = SPivot
- params['R'] = SPivot
- params['SMarginal'] = SMarginal
- params['MMarginal'] = MMarginal
- params['polybasisMarginal'] = "MONOMIAL" + radialM
- params['radialDirectionalWeightsMarginal'] = radialWeightM
- params['matchingWeight'] = matchingWeight
- params['cutOffTolerance'] = cutOffTolerance
- params["cutOffType"] = cutOffType
- method = RBP
+ raise
+# params['S'] = SPivot
+# params['R'] = SPivot
+# params['SMarginal'] = SMarginal
+# params['MMarginal'] = MMarginal
+# params['polybasisMarginal'] = "MONOMIAL" + radialM
+# params['radialDirectionalWeightsMarginal'] = radialWeightM
+# params['matchingWeight'] = matchingWeight
+# params['cutOffTolerance'] = cutOffTolerance
+# params["cutOffType"] = cutOffType
+# method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
- params['S'] = STensorized
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": 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', what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF', what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err', what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res', 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)))
approx.verbosity = 5
-from plot_zero_set import plotZeroSet2
-muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
- 200, [2., 1.])
+try:
+ from plot_zero_set import plotZeroSet2
+ muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
+ 200, [2., 1.])
+except:
+ pass
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 25,
[2., 1.], relative = True,
nobeta = True)
-print(np.sort(approx.getPoles([None, .5]) ** 2.))
+try:
+ print(np.sort(approx.getPoles([None, .5]) ** 2.))
+except:
+ pass
diff --git a/examples/2d/pod/matrix_passive_pod.py b/examples/2d/pod/matrix_passive_pod.py
index 8c3b8cb..1789f6f 100644
--- a/examples/2d/pod/matrix_passive_pod.py
+++ b/examples/2d/pod/matrix_passive_pod.py
@@ -1,192 +1,193 @@
import numpy as np
from rrompy.hfengines.linear_problem.tridimensional import \
MatrixDynamicalPassive as MDP
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
-from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
+#from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 10
size = 3
show_sample = True
show_norm = True
Delta = 0
MN = 5
R = (MN + 2) * (MN + 1) // 2
STensorized = (MN + 1) ** 2
PODTol = 1e-6
SPivot = MN + 1
SMarginal = 3
MMarginal = SMarginal - 1
samples = "centered"
samples = "standard"
samples = "pivoted"
algo = "rational"
algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples in ["standard", "pivoted"]:
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 10.
radialWeight = [rW0]
if samples == "pivoted":
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
matchingWeight = 1.
cutOffTolerance = 5.
cutOffType = "POTENTIAL"
if size == 1:
mu0 = [2.7e2, 20]
mutar = [3e2, 25]
murange = [[20., 10], [5.2e2, 30]]
elif size == 2:
mu0 = [2.7e2, 60]
mutar = [3e2, 75]
murange = [[20., 10], [5.2e2, 110]]
elif size == 3:
mu0 = [2.7e2, 160]
mutar = [3e2, 105]
murange = [[20., 10], [5.2e2, 310]]
assert Delta <= 0
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[aEff*murange[0][0] + bEff*murange[1][0],
aEff*murange[0][1] + bEff*murange[1][1]],
[aEff*murange[1][0] + bEff*murange[0][0],
aEff*murange[1][1] + bEff*murange[0][1]]]
n = 100
b = 10
solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb)
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight * 2
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "pivoted":
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
params['matchingWeight'] = matchingWeight
params['cutOffTolerance'] = cutOffTolerance
params["cutOffType"] = cutOffType
method = RIP
else: #if algo == "RB":
params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':R,
'POD':True, 'PODTolerance':PODTol}
if samples == "standard":
method = RB
elif samples == "centered":
method = RB
elif samples == "pivoted":
- params['R'] = SPivot
- params['S'] = SPivot
- params['SMarginal'] = SMarginal
- params['MMarginal'] = MMarginal
- params['polybasisMarginal'] = "MONOMIAL" + radialM
- params['radialDirectionalWeightsMarginal'] = radialWeightM
- params['matchingWeight'] = matchingWeight
- params['cutOffTolerance'] = cutOffTolerance
- params["cutOffType"] = cutOffType
- method = RBP
+ raise
+# params['R'] = SPivot
+# params['S'] = SPivot
+# params['SMarginal'] = SMarginal
+# params['MMarginal'] = MMarginal
+# params['polybasisMarginal'] = "MONOMIAL" + radialM
+# params['radialDirectionalWeightsMarginal'] = radialWeightM
+# params['matchingWeight'] = matchingWeight
+# params['cutOffTolerance'] = cutOffTolerance
+# params["cutOffType"] = cutOffType
+# method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV")
# params['sampler'] = QS(murange, "GAUSSLEGENDRE")
# params['sampler'] = QS(murange, "UNIFORM")
params["S"] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV")
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON")
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0])
elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
L = mutar[1]
approx.plotApprox(mutar, name = 'u_app', what = "REAL")
approx.plotHF(mutar, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, name = 'err', what = "REAL")
# approx.plotRes(mutar, name = 'res', 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)))
try:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [1., 1], polesImTol = 2.)
except: pass
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[1., 1.], relative = False,
nobeta = True)
print(1.j * approx.getPoles([None, 50.]))
diff --git a/examples/2d/pod/plot_inf_set.py b/examples/2d/pod/plot_inf_set.py
index 00aa698..a19285e 100644
--- a/examples/2d/pod/plot_inf_set.py
+++ b/examples/2d/pod/plot_inf_set.py
@@ -1,362 +1,362 @@
import warnings
import numpy as np
from matplotlib import pyplot as plt
def plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
exp = 2., normalizeDen = False):
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)
if not np.isnan(beta[0]):
eta = R / beta / E
betamin, betamax = np.min(beta), np.max(beta)
etamin, etamax = np.min(eta), np.max(eta)
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(np.real(mu2x), [ZTmin] * len(mu2x), 'kx')
for j, x in enumerate(np.real(mu2x)):
plt.annotate("{}".format(j), (x, ZTmin))
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(np.real(mu2x), [Rmax] * len(mu2x), 'kx')
for j, x in enumerate(np.real(mu2x)):
plt.annotate("{}".format(j), (x, Rmax))
plt.plot([murangeExp[0][0]] * 2, [Rmin, Rmax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Rmin, Rmax], 'm:')
plt.xlim(mu1[0], mu1[-1])
if normalizeDen:
plt.title("|Q(mu)res(mu)|")
else:
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(np.real(mu2x), [Emax] * len(mu2x), 'kx')
for j, x in enumerate(np.real(mu2x)):
plt.annotate("{}".format(j), (x, Emax))
plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
plt.xlim(mu1[0], mu1[-1])
if normalizeDen:
plt.title("|Q(mu)err(mu)|")
else:
plt.title("|err(mu)|")
plt.grid()
plt.show()
if not np.isnan(beta[0]):
plt.figure(figsize = (15, 7))
plt.jet()
plt.plot(mu1, beta)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
plt.plot(np.real(mu2x), [betamax] * len(mu2x), 'kx')
for j, x in enumerate(np.real(mu2x)):
plt.annotate("{}".format(j), (x, betamax))
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()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, R / beta)
plt.semilogy(mu1, E, '--')
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
plt.plot(np.real(mu2x), [Emax] * len(mu2x), 'kx')
for j, x in enumerate(np.real(mu2x)):
plt.annotate("{}".format(j), (x, Emax))
plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
plt.xlim(mu1[0], mu1[-1])
if normalizeDen:
plt.title("|Q(mu)res(mu)/beta(mu)|")
else:
plt.title("|res(mu)/beta(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, eta)
for l_ in approx.trainedModel.getPoles():
plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
plt.plot(np.real(mu2x), [etamax] * len(mu2x), 'kx')
for j, x in enumerate(np.real(mu2x)):
plt.annotate("{}".format(j), (x, etamax))
plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("eta(mu)")
plt.grid()
plt.show()
def plotInfSet1(murange, murangeEff, approx, mu0, nSamples = 20, exp = 2.,
relative = True, normalizeDen = False, nobeta = False):
mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
nSamples)
mus = np.power(mu1, 1. / exp)
Z = approx.normHF(mus)
T = approx.normApprox(mus)
R = approx.normRes(mus)
E = approx.normErr(mus)
if relative:
F = approx.normRHS(mus)
R /= F
E /= Z
if normalizeDen:
Qvals = np.abs(approx.trainedModel.getQVal(mus))
R *= Qvals
E *= Qvals
if nobeta:
beta = np.empty(len(mus))
beta[:] = np.nan
else:
- beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
+ beta = approx.HFEngine.stabilityFactor(mus, approx.getHF(mus))
plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
exp, normalizeDen)
return mus, Z, T, R, E, beta
def plotInfSet2FromData(mus, Ze, Te, Re, Ee, beta, murange, approx, mu0,
exps = [2., 2.], clip = -1, normalizeDen = False):
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.log10(Ze)
T = np.log10(Te)
R = np.log10(Re)
E = np.log10(Ee)
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)
if not np.isnan(beta[0, 0]):
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,
np.ComplexWarning))
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(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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("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(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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|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(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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:')
if normalizeDen:
plt.title("log10|Q(mu)res(mu)|")
else:
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(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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:')
if normalizeDen:
plt.title("log10|Q(mu)err(mu)|")
else:
plt.title("log10|err(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
if not np.isnan(beta[0, 0]):
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, beta,
levels = np.linspace(betamin, betamax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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("beta(mu)")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
levels = np.linspace(Emin, Emax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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:')
if normalizeDen:
plt.title("log10|Q(mu)res(mu)/beta(mu)|")
else:
plt.title("log10|res(mu)/beta(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 50)
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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|eta(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
def plotInfSet2(murange, murangeEff, approx, mu0, nSamples = 200,
exps = [2., 2.], clip = -1, relative = True,
normalizeDen = False, nobeta = False):
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]
Ze = approx.normHF(mus).reshape((nSamples, nSamples))
Te = approx.normApprox(mus).reshape((nSamples, nSamples))
Re = approx.normRes(mus).reshape((nSamples, nSamples))
Ee = approx.normErr(mus).reshape((nSamples, nSamples))
if relative:
Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
Re /= Fe
Ee /= Ze
if normalizeDen:
Qvals = np.abs(approx.trainedModel.getQVal(mus).reshape(
(nSamples, nSamples)))
Re *= Qvals
Ee *= Qvals
if nobeta:
betae = np.empty((nSamples, nSamples))
betae[:, :] = np.nan
else:
- betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)\
+ betae = approx.HFEngine.stabilityFactor(mus, approx.getHF(mus))\
.reshape((nSamples, nSamples))
plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange,
approx, mu0, exps, clip, normalizeDen)
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 8e82cf0..e62acd0 100644
--- a/examples/2d/pod/plot_zero_set.py
+++ b/examples/2d/pod/plot_zero_set.py
@@ -1,99 +1,103 @@
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)
+ try:
+ Z = approx.trainedModel.getQVal(mus)
+ Zabs = np.abs(Z)
+ Zmin, Zmax = np.min(Zabs), np.max(Zabs)
+ except:
+ Z = None
+ Zmin, Zmax = 0., 1.
plt.figure(figsize = (15, 7))
plt.jet()
- plt.semilogy(mu1, Zabs)
+ if Z is not None:
+ 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')
for j, x in enumerate(mu2x):
plt.annotate("{}".format(j), (x, Zmax))
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, polesImTol : float = 1e-5):
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]
poles1, poles2 = [], []
for m2 in mu2s:
pls = approx.getPoles([None, m2]) ** exps[0]
pls[np.imag(pls) > polesImTol] = np.nan
pls = np.real(pls)
poles1 += list(pls)
poles2 += [m2] * len(pls)
- Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
- 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
try:
- pass
+ Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape)
+ 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
except:
Z = None
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = (UserWarning,
np.ComplexWarning))
plt.figure(figsize = (15, 7))
plt.jet()
if Z is not None:
p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap,
levels = np.linspace(Zabsmin, Zabsmax, 50))
if clip > 0:
plt.contour(Mu1, Mu2, Zabs, [Zabsmin])
plt.plot(poles1, poles2, 'k.')
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
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:')
if Z is not None:
plt.colorbar(p)
plt.xlim(murangeExp[0][0], murangeExp[1][0])
plt.ylim(murangeExp[0][1], murangeExp[1][1])
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 fd4cbf7..2d11ffa 100644
--- a/examples/2d/pod/square_pod.py
+++ b/examples/2d/pod/square_pod.py
@@ -1,187 +1,192 @@
import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
HelmholtzSquareDomainProblemEngine as HSDPE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
-from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
+#from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
-size = 7
+size = 6
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
STensorized = (MN + 1) ** 2
PODTol = 1e-6
SPivot = MN + 1
-SMarginal = 3
+SMarginal = 4
MMarginal = SMarginal - 1
-matchingWeight = 10.
samples = "centered"
samples = "standard"
samples = "pivoted"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
#samplingM = "quadrature_total"
#samplingM = "random"
if samples == "pivoted":
radialM = ""
radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0]
+ cutOffTolerance = 1. * np.inf
+ cutOffType = "MAGNITUDE"
+ matchingWeight = 10.
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.1
+aEff = 1.#1
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
rescalingExp = [2., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV"
# params['polybasis'] = "LEGENDRE"
# params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "pivoted":
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV"
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
+ params['cutOffTolerance'] = cutOffTolerance
+ params["cutOffType"] = cutOffType
method = RIP
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
- if samples == "standard":
- method = RB
- elif samples == "centered":
+ if samples in ["standard", "centered"]:
method = RB
elif samples == "pivoted":
- params['S'] = SPivot
- params['SMarginal'] = SMarginal
- params['MMarginal'] = MMarginal
- params['polybasisMarginal'] = "MONOMIAL" + radialM
- params['matchingWeight'] = matchingWeight
- params['radialDirectionalWeightsMarginal'] = radialWeightM
- method = RBP
+ raise
+# params['S'] = SPivot
+# params['SMarginal'] = SMarginal
+# params['MMarginal'] = MMarginal
+# params['polybasisMarginal'] = "MONOMIAL" + radialM
+# params['matchingWeight'] = matchingWeight
+# params['radialDirectionalWeightsMarginal'] = radialWeightM
+# params['cutOffTolerance'] = cutOffTolerance
+# params["cutOffType"] = cutOffType
+# method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1], murange[1][1]], "UNIFORM")
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1], murange[1][1]], "CHEBYSHEV")
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1], murange[1][1]], "HALTON")
if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
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
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 20, [2., 1.], clip = clip,
- relative = False, nobeta = True)
+ relative = True, nobeta = True)
diff --git a/examples/3d/base/fracture3.py b/examples/3d/base/fracture3.py
index b987ab8..13b6a32 100644
--- a/examples/3d/base/fracture3.py
+++ b/examples/3d/base/fracture3.py
@@ -1,35 +1,35 @@
from rrompy.hfengines.linear_problem.tridimensional import \
MembraneFractureEngine3 as MFE
from fracture3_warping import fracture3_warping
verb = 100
mu0 = [45. ** .5, .7, .075]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta,
n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu0)[0]
#solver.plotmesh(figsize = (7.5, 4.5))
#solver.plot(u0, what = 'REAL', figsize = (8, 5))
print(solver.norm(uh))
#solver.plot(uh, what = 'REAL', figsize = (8, 5))
-#solver.plot(solver.residual(uh, mu0)[0], name = 'res',
+#solver.plot(solver.residual(mu0, uh)[0], name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], filename = 'out', folder = True,
# forceNewFile = False)
warps = fracture3_warping(solver.V.mesh(), L, mu0[1], delta, mu0[2])
#solver.plotmesh(warps, figsize = (7.5, 4.5))
#solver.plot(u0, warps, what = 'REAL', figsize = (8, 5))
solver.plot(uh, warps, what = 'REAL', figsize = (8, 5))
-#solver.plot(solver.residual(uh, mu0)[0], warps, name = 'res',
+#solver.plot(solver.residual(mu0, uh)[0], warps, name = 'res',
# what = 'REAL', figsize = (8, 5))
#solver.outParaviewTimeDomain(uh, mu0[0], warps, filename = 'outW',
# folder = True, forceNewFile = False)
diff --git a/examples/3d/base/scattering1.py b/examples/3d/base/scattering1.py
index 0f304b0..f9a0c2d 100644
--- a/examples/3d/base/scattering1.py
+++ b/examples/3d/base/scattering1.py
@@ -1,29 +1,29 @@
from rrompy.hfengines.linear_problem.tridimensional import Scattering1d as S1D
import fenics as fen
import numpy as np
verb = 100
mu0 = [4., np.pi, 0.]
mu = [4.5, np.pi, 1.]
n = 100
solver = S1D(mu0 = mu0, n = n, verbosity = verb)
u0 = solver.liftDirichletData()
uh = solver.solve(mu)[0]
#solver.plotmesh(figsize = (12, 2))
#solver.plot(u0, what = ['REAL', 'IMAG'], figsize = (12, 4))
print(solver.norm(uh))
#solver.plot(uh, what = ['REAL', 'IMAG'], figsize = (12, 4))
-#solver.plot(solver.residual(uh, mu0)[0], name = 'res',
+#solver.plot(solver.residual(mu0, uh)[0], name = 'res',
# what = ['REAL', 'IMAG'], figsize = (12, 4))
x = fen.SpatialCoordinate(solver.V.mesh())
warps = [x * mu[1] / solver._L - x, x * solver._L / mu[1] - x]
#solver.plotmesh(warps, figsize = (12, 2))
#solver.plot(u0, warps, what = ['REAL', 'IMAG'], figsize = (12, 4))
solver.plot(uh, warps, what = ['REAL', 'IMAG'], figsize = (12, 4))
-#solver.plot(solver.residual(uh, mu0)[0], warps, name = 'res',
+#solver.plot(solver.residual(mu0, uh)[0], warps, name = 'res',
# what = ['REAL', 'IMAG'], figsize = (12, 4))
diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py
index 7ec301c..639f574 100644
--- a/examples/3d/pod/fracture3_pod.py
+++ b/examples/3d/pod/fracture3_pod.py
@@ -1,242 +1,243 @@
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
from rrompy.hfengines.linear_problem.tridimensional import \
MembraneFractureEngine3 as MFE
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted as RIP
-from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
+#from rrompy.reduction_methods.pivoted import ReducedBasisPivoted as RBP
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 70
size = 4
show_sample = False
show_norm = False
clip = -1
#clip = .4
#clip = .6
Delta = 0
MN = 5
R = (MN + 3) * (MN + 2) * (MN + 1) // 6
STensorized = (MN + 1) ** 3
PODTol = 1e-8
SPivot = MN + 1
SMarg1d = 5
MMarginal = SMarg1d - 1
SMarginal = (MMarginal + 2) * (MMarginal + 1) // 2
SMarginalTensorized = SMarg1d ** 2
matchingWeight = 10.
samples = "centered"
samples = "standard"
#samples = "pivoted"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
#sampling = "quadrature_total"
#sampling = "random"
samplingM = "quadrature"
samplingM = "quadrature_total"
#samplingM = "random"
polydegreetype = "TOTAL"
#polydegreetype = "FULL"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 3
if samples == "pivoted":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0]
radialM = ""
# radialM = "_gaussian"
# radialM = "_thinplate"
# radialM = "_multiquadric"
rMW0 = 2.
radialWeightM = [rMW0] * 2
assert Delta <= 0
if size == 1:
mu0 = [47.5 ** .5, .4, .05]
mutar = [50 ** .5, .45, .07]
murange = [[40 ** .5, .3, .025], [55 ** .5, .5, .075]]
if size == 2:
mu0 = [50 ** .5, .3, .02]
mutar = [55 ** .5, .4, .03]
murange = [[40 ** .5, .1, .005], [60 ** .5, .5, .035]]
if size == 3:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .4, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
if size == 4:
mu0 = [45 ** .5, .4, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .5, .06]]
if size == 5:
mu0 = [45 ** .5, .5, .05]
mutar = [47 ** .5, .45, .045]
murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
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[0][2] + bEff*murange[1][2]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1],
aEff*murange[1][2] + bEff*murange[0][2]]]
H = 1.
L = .75
delta = .05
n = 50
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1., 1.]
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True,
'polydegreetype':polydegreetype}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
method = RI
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
elif samples == "pivoted":
params['S'] = SPivot
params['SMarginal'] = SMarginal
params['MMarginal'] = MMarginal
params['polybasisPivot'] = "CHEBYSHEV" + radial
params['polybasisMarginal'] = "MONOMIAL" + radialM
params['matchingWeight'] = matchingWeight
params['radialDirectionalWeightsPivot'] = radialWeight
params['radialDirectionalWeightsMarginal'] = radialWeightM
method = RIP
else: #if algo == "RB":
params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
'S':R, 'POD':True, 'PODTolerance':PODTol}
if samples in ["standard", "centered"]:
method = RB
elif samples == "pivoted":
- params['S'] = SPivot
- params['SMarginal'] = SMarginal
- params['MMarginal'] = MMarginal
- params['polybasisMarginal'] = "MONOMIAL" + radialM
- params['matchingWeight'] = matchingWeight
- params['radialDirectionalWeightsMarginal'] = radialWeightM
- method = RBP
+ raise
+# params['S'] = SPivot
+# params['SMarginal'] = SMarginal
+# params['MMarginal'] = MMarginal
+# params['polybasisMarginal'] = "MONOMIAL" + radialM
+# params['matchingWeight'] = matchingWeight
+# params['radialDirectionalWeightsMarginal'] = radialWeightM
+# method = RBP
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE",scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
params['S'] = STensorized
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
elif samples == "pivoted":
if sampling == "quadrature":
params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "CHEBYSHEV")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "GAUSSLEGENDRE")
# params['samplerPivot'] = QS([murange[0][0], murange[1][0]], "UNIFORM")
elif sampling == "quadrature_total":
params['samplerPivot'] = QST([murange[0][0], murange[1][0]], "CHEBYSHEV")
else: # if sampling == "random":
params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON")
if samplingM == "quadrature":
params['samplerMarginal'] = QS([murange[0][1:], murange[1][1:]], "UNIFORM")
params['SMarginal'] = SMarginalTensorized
elif samplingM == "quadrature_total":
params['samplerMarginal'] = QST([murange[0][1:], murange[1][1:]], "CHEBYSHEV")
params['polybasisMarginal'] = "CHEBYSHEV" + radialM
else: # if samplingM == "random":
params['samplerMarginal'] = RS([murange[0][1:], murange[1][1:]], "HALTON")
if samples != "pivoted":
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
else:
approx = method(solver, mu0 = mu0, directionPivot = [0],
approxParameters = params, verbosity = verb)
if samples != "centered": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
from fracture3_warping import fracture3_warping
warps = fracture3_warping(solver.V.mesh(), L, mutar[1], delta, mutar[2])
approx.plotApprox(mutar, warps, name = 'u_app', what = "REAL")
approx.plotHF(mutar, warps, name = 'u_HF', what = "REAL")
approx.plotErr(mutar, warps, name = 'err', what = "REAL")
# approx.plotRes(mutar, warps, name = 'res', 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)))
fig = plt.figure(figsize = (8, 6))
ax = Axes3D(fig)
ax.scatter(np.real(approx.mus(0) ** 2.), np.real(approx.mus(1)),
np.real(approx.mus(2)), '.')
plt.show()
plt.close()
approx.verbosity = 0
approx.trainedModel.verbosity = 0
from plot_zero_set_3 import plotZeroSet3
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, None, mu0[2]], [2., 1., 1.],
# clip = clip)
#muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
# [None, mu0[1], None], [2., 1., 1.],
# clip = clip)
muZeroScatter = plotZeroSet3(murange, murangeEff, approx, mu0, 50,
[None, None, None], [2., 1., 1.], clip = clip)
if show_norm:
solver._solveBatchSize = 25
from plot_inf_set_3 import plotInfSet3
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, None, mu0[2]], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
murange, murangeEff, approx, mu0, 25,
[None, mu0[1], None], [2., 1., 1.],
clip = clip, relative = False,
normalizeDen = True)
print(approx.getPoles([None, .5, .05]) ** 2.)
diff --git a/examples/3d/pod/plot_inf_set_3.py b/examples/3d/pod/plot_inf_set_3.py
index c823ed3..900e447 100644
--- a/examples/3d/pod/plot_inf_set_3.py
+++ b/examples/3d/pod/plot_inf_set_3.py
@@ -1,360 +1,360 @@
import warnings
from copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
def plotInfSet31FromData(mNone, mus, Z, T, R, E, beta, murange, approx, mu0,
exp = 2., normalizeDen = False):
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exp
else:
mu2x = mu0[0] ** exp
murangeExp = [[murange[0][mNone[0]] ** exp],
[murange[1][mNone[0]] ** exp]]
mu1 = np.real(np.power([m[mNone[0]] for m in 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)
if not np.isnan(beta[0]):
eta = R / beta / E
betamin, betamax = np.min(beta), np.max(beta)
etamin, etamax = np.min(eta), np.max(eta)
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, Z)
plt.semilogy(mu1, T, '--')
mTMP = mus[0]
mTMP[mNone[0]] = None
for l_ in approx.trainedModel.getPoles(mTMP):
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(mTMP):
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])
if normalizeDen:
plt.title("|Q(mu)res(mu)|")
else:
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(mTMP):
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])
if normalizeDen:
plt.title("|Q(mu)err(mu)|")
else:
plt.title("|err(mu)|")
plt.grid()
plt.show()
if not np.isnan(beta[0]):
plt.figure(figsize = (15, 7))
plt.jet()
plt.plot(mu1, beta)
for l_ in approx.trainedModel.getPoles(mTMP):
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()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, R / beta)
plt.semilogy(mu1, E, '--')
for l_ in approx.trainedModel.getPoles(mTMP):
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])
if normalizeDen:
plt.title("|Q(mu)res(mu)/beta(mu)|")
else:
plt.title("|res(mu)/beta(mu)|")
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
plt.semilogy(mu1, eta)
for l_ in approx.trainedModel.getPoles(mTMP):
plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
plt.plot(mu2x, [etamax] * len(mu2x), 'kx')
plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
plt.xlim(mu1[0], mu1[-1])
plt.title("eta(mu)")
plt.grid()
plt.show()
def plotInfSet32FromData(mNone, mus, Ze, Te, Re, Ee, beta, murange, approx,
mu0, exps = [2., 1.], clip = -1,
normalizeDen = False):
if hasattr(approx, "mus"):
mu2x = approx.mus(mNone[0]) ** exps[0]
mu2y = approx.mus(mNone[1]) ** exps[1]
else:
mu2x, mu2y = mu0[mNone[0]] ** exps[0], mu0[mNone[1]] ** exps[1]
murangeExp = [[murange[0][mNone[0]] ** exps[0],
murange[0][mNone[1]] ** exps[1]],
[murange[1][mNone[0]] ** exps[0],
murange[1][mNone[1]] ** exps[1]]]
mu1s = np.unique([m[mNone[0]] for m in mus])
mu2s = np.unique([m[mNone[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.log10(Ze)
T = np.log10(Te)
R = np.log10(Re)
E = np.log10(Ee)
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)
if not np.isnan(beta[0, 0]):
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,
np.ComplexWarning))
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("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("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:')
if normalizeDen:
plt.title("log10|Q(mu)res(mu)|")
else:
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:')
if normalizeDen:
plt.title("log10|Q(mu)err(mu)|")
else:
plt.title("log10|err(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
if not np.isnan(beta[0, 0]):
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("beta(mu)")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
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:')
if normalizeDen:
plt.title("log10|Q(mu)res(mu)/beta(mu)|")
else:
plt.title("log10|res(mu)/beta(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 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|eta(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
def plotInfSet3(murange, murangeEff, approx, mu0, nSamples = 200,
marginal = [None, None, .05], exps = [2., 1., 1.], clip = -1,
relative = True, normalizeDen = False, nobeta = False):
mNone = [i for i, m in enumerate(marginal) if m is None]
nNone = len(mNone)
if nNone not in [1, 2]:
raise
if nNone == 1:
exp = exps[mNone[0]]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exp,
murangeEff[1][mNone[0]] ** exp, nSamples)
mu1s = np.power(mu1, 1. / exp)
mus = []
mBase = copy(marginal)
for m1 in mu1s:
mBase[mNone[0]] = m1
mus += [copy(mBase)]
Z = approx.normHF(mus)
T = approx.normApprox(mus)
R = approx.normRes(mus)
E = approx.normErr(mus)
if relative:
F = approx.normRHS(mus)
R /= F
E /= Z
if normalizeDen:
Qvals = np.abs(approx.trainedModel.getQVal(mus))
R *= Qvals
E *= Qvals
if nobeta:
beta = np.empty(len(mus))
beta[:] = np.nan
else:
- beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
+ beta = approx.HFEngine.stabilityFactor(mus, approx.getHF(mus))
plotInfSet31FromData(mNone, mus, Z, T, R, E, beta, murange, approx,
mu0, exp, normalizeDen)
return mus, Z, T, R, E, beta
if nNone == 2:
exps = [exps[m] for m in mNone]
mu1 = np.linspace(murangeEff[0][mNone[0]] ** exps[0],
murangeEff[1][mNone[0]] ** exps[0], nSamples)
mu2 = np.linspace(murangeEff[0][mNone[1]] ** exps[1],
murangeEff[1][mNone[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 = []
mBase = copy(marginal)
for m2 in mu2s:
mBase[mNone[1]] = m2
for m1 in mu1s:
mBase[mNone[0]] = m1
mus += [copy(mBase)]
Ze = approx.normHF(mus).reshape((nSamples, nSamples))
Te = approx.normApprox(mus).reshape((nSamples, nSamples))
Re = approx.normRes(mus).reshape((nSamples, nSamples))
Ee = approx.normErr(mus).reshape((nSamples, nSamples))
if relative:
Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
Re /= Fe
Ee /= Ze
if normalizeDen:
Qvals = np.abs(approx.trainedModel.getQVal(mus).reshape(
(nSamples, nSamples)))
Re *= Qvals
Ee *= Qvals
if nobeta:
betae = np.empty((nSamples, nSamples))
betae[:, :] = np.nan
else:
- betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)\
+ betae = approx.HFEngine.stabilityFactor(mus, approx.getHF(mus))\
.reshape((nSamples, nSamples))
plotInfSet32FromData(mNone, mus, Ze, Te, Re, Ee, betae, murange,
approx, mu0, exps, clip, normalizeDen)
return mus, Ze, Te, Re, Ee, betae
diff --git a/examples/all_forcing/all_forcing_engine.py b/examples/all_forcing/all_forcing_engine.py
index 94980c1..72fa63d 100644
--- a/examples/all_forcing/all_forcing_engine.py
+++ b/examples/all_forcing/all_forcing_engine.py
@@ -1,42 +1,43 @@
import numpy as np
import fenics as fen
from rrompy.hfengines.linear_problem import LaplaceBaseProblemEngine as LBPE
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.solver.fenics import fenics2Vector
class AllForcingEngine(LBPE):
def __init__(self, mu0:float, n : int = 30,
degree_threshold : int = np.inf, verbosity : int = 10,
timestamp : bool = True):
super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
homogeneized = False, verbosity = verbosity,
timestamp = timestamp)
mesh = fen.RectangleMesh(fen.Point(-5., -5.), fen.Point(5., 5.), n, n)
self.V = fen.FunctionSpace(mesh, "P", 1)
self.nAs, self.nbs = 1, 4
x, y = fen.SpatialCoordinate(mesh)[:]
scaling = (2. * np.pi) ** -1.
r2 = x ** 2. + y ** 2.
self.forcingCoeffs = [
scaling * fen.exp(- (r2 + 1. - 2. * x + 1. - 2. * y) / 2. / 4.) / 2.,
scaling * fen.exp(- (r2 + 1. + 2. * x + 1. + 2. * y) / 2. / 16.) / 4.,
- scaling * fen.exp(- (r2 + 1. + 2. * x + 1. - 2. * y) / 2. / 9.) / 30.,
scaling * fen.exp(- (r2 + 1. - 2. * x + 1. + 2. * y) / 2. / 25.) / 120.]
def b(self, mu = [], der = 0):
- mu = self.checkParameter(mu)
- if not hasattr(der, "__len__"): der = [der] * self.npar
- derI = der[0]
- for j in range(derI, self.nbs):
+ derI = hashD(der) if hasattr(der, "__len__") else der
+ if self.thbs[0] is None:
+ self.thbs = [None] * (self.nbs + self.homogeneized * self.nAs)
+ self._setbHomogeneized()
+ for j in range(self.nbs):
if self.bs[j] is None:
vbMng(self, "INIT", "Assembling forcing term b{}.".format(j),
20)
parsRe = self.iterReduceQuadratureDegree([(
self.forcingCoeffs[j],
"forcingCoefficient")])
u0Re = self.DirichletDatum[0]
L0Re = fen.dot(self.forcingCoeffs[j], self.v) * fen.dx
DBCR = fen.DirichletBC(self.V, u0Re, self.DirichletBoundary)
self.bs[j] = fenics2Vector(L0Re, parsRe, DBCR, 1)
vbMng(self, "DEL", "Done assembling forcing term.", 20)
- return self._assembleb(mu, der, derI)
+ return self._assembleb(mu, derI)
diff --git a/examples/base/helmholtz_resonator.py b/examples/base/helmholtz_resonator.py
index d5d1d9a..9a9e441 100644
--- a/examples/base/helmholtz_resonator.py
+++ b/examples/base/helmholtz_resonator.py
@@ -1,40 +1,44 @@
from matplotlib import pyplot as plt
import fenics as fen
import mshr
import ufl
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
p = plt.jet()
n = 50
boundary = mshr.Polygon([fen.Point(0, 0), fen.Point(6, 0), fen.Point(6, 1),
fen.Point(1.3, 1), fen.Point(1.3, 1.2),
fen.Point(1.65, 1.2), fen.Point(1.65, 2.2),
fen.Point(.65, 2.2), fen.Point(.65, 1.2),
fen.Point(1, 1.2), fen.Point(1, 1), fen.Point(0, 1)])
mesh = mshr.generate_mesh(boundary, n)
-for k in range(10):
+for k in range(1):
kappa = .25 + .05 * k
ZR = 10.
x, y = fen.SpatialCoordinate(mesh)[:]
solver = HPE(kappa)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.RobinBoundary = (lambda x, on_b: on_b and (fen.near(x[0], 6.)
- or fen.near(x[1], 2.25)))
+ or fen.near(x[1], 2.2)))
solver.NeumannBoundary = "REST"
solver.signR = + 1.
solver.NeumannDatum = [fen.Constant(0.),
ufl.conditional(ufl.And(ufl.gt(y, 0.), ufl.lt(y, 1.)),
fen.Constant(kappa), fen.Constant(0.))]
solver.RobinDatumH = [fen.Constant(0.),
ufl.conditional(ufl.gt(y, 1.25),
fen.Constant(kappa / ZR),
fen.Constant(kappa))]
uh = solver.solve(kappa)[0]
solver.plot(uh, name = "k={}".format(kappa))
print(solver.norm(uh))
+from rrompy.hfengines.base import MarginalProxyEngine
+
+ms = MarginalProxyEngine(solver, [None])
+print(ms.A(1))
\ No newline at end of file
diff --git a/examples/base/matrix_solver.py b/examples/base/matrix_solver.py
index 5db946c..536ef9f 100644
--- a/examples/base/matrix_solver.py
+++ b/examples/base/matrix_solver.py
@@ -1,30 +1,30 @@
import numpy as np
import scipy.sparse as sp
from rrompy.hfengines.base import MatrixEngineBase as MEB
test = 2
N = 100
verb = 0
solver = MEB(verbosity = verb)
solver.npar = 1
solver.nAs = 2
mu = 5.001
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)]
solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
uh = solver.solve(mu)[0]
print(solver.norm(uh))
solver.plot(uh)
-solver.plot(solver.residual(uh, mu)[0], 'res')
+solver.plot(solver.residual(mu, uh)[0], 'res')
diff --git a/examples/base/solver.py b/examples/base/solver.py
index 1217af3..a0be1bd 100644
--- a/examples/base/solver.py
+++ b/examples/base/solver.py
@@ -1,75 +1,75 @@
import numpy as np
from rrompy.hfengines.linear_problem import \
HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.linear_problem import \
HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.hfengines.linear_problem import \
HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.hfengines.linear_problem import \
HelmholtzCavityScatteringProblemEngine as HCSPE
-testNo = 4
+testNo = 2
verb = 0
if testNo == 1:
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20,
verbosity = verb)
mu = 12.**.5
solver.setSolver("BICG", {"tol" : 1e-15})
uh = solver.solve(mu)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
- solver.plot(solver.residual(uh, mu)[0], 'res')
+ solver.plot(solver.residual(mu, uh)[0], name = 'res')
###########
elif testNo in [2, -2]:
solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4.,
- n = 50, verbosity = verb, homogeneized = (testNo < 0))
+ n = 50, verbosity = verb)
mu = 4.
uref = solver.liftDirichletData()
if testNo > 0:
uh = solver.solve(mu)[0]
utot = uh - uref
else:
utot = solver.solve(mu)[0]
uh = utot + uref
print(solver.norm(uh))
print(solver.norm(uref))
solver.plot(uh)
solver.plot(uref, name = 'u_Dir')
solver.plot(utot, name = 'u_tot')
- solver.plot(solver.residual(uh, mu)[0], 'res')
- solver.plot(solver.residual(utot, mu)[0], 'res_tot')
+ solver.plot(solver.residual(mu, uh)[0], name = 'res')
+ solver.plot(solver.residual(mu, utot)[0], name = 'res_tot')
###########
elif testNo in [3, -3]:
solver = HBSPE(R = 5, kappa = 12**.5, theta = - np.pi * 60 / 180, n = 30,
verbosity = verb, homogeneized = (testNo < 0))
mu = 12**.5
uref = solver.liftDirichletData()
if testNo > 0:
uh = solver.solve(mu)[0]
utot = uh - uref
else:
utot = solver.solve(mu)[0]
uh = utot + uref
solver.plotmesh()
print(solver.norm(uh))
print(solver.norm(utot))
solver.plot(uh)
solver.plot(utot, name = 'u_tot')
- solver.plot(solver.residual(uh, mu)[0], 'res')
- solver.plot(solver.residual(utot, mu)[0], 'res_tot')
+ solver.plot(solver.residual(mu, uh)[0], 'res')
+ solver.plot(solver.residual(mu, utot)[0], 'res_tot')
###########
elif testNo == 4:
solver = HCSPE(kappa = 5, n = 30, verbosity = verb)
mu = 10
uh = solver.solve(mu)[0]
solver.plotmesh()
print(solver.norm(uh))
solver.plot(uh)
- solver.plot(solver.residual(uh, mu)[0], 'res')
+ solver.plot(solver.residual(mu, uh)[0], 'res')
diff --git a/examples/from_papers/greedy_internalBox.py b/examples/from_papers/greedy_internalBox.py
index 336f8d5..adc9ea4 100644
--- a/examples/from_papers/greedy_internalBox.py
+++ b/examples/from_papers/greedy_internalBox.py
@@ -1,117 +1,117 @@
import numpy as np
import fenics as fen
from rrompy.hfengines.linear_problem import HelmholtzProblemEngine as HPE
from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RI
from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RB
from rrompy.solver.fenics import L2NormMatrix
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
dim = 2
verb = 5
timed = False
algo = "RI"
#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 = {'sampler':QS([kl, kr], "UNIFORM", 2.), 'nTestPoints':500,
- 'greedyTol':1e2, 'S':2, 'polybasis':polyBasis,
+ 'greedyTol':1e-2, 'S':2, 'polybasis':polyBasis,
# 'errorEstimatorKind':'AFFINE'}
'errorEstimatorKind':'DISCREPANCY'}
# 'errorEstimatorKind':'INTERPOLATORY'}
# 'errorEstimatorKind':'LOCALL2'}
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(np.real(k0), verbosity = verb)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.refractionIndex = fen.Constant(1. / 54.6)
solver.forcingTerm = f
solver.NeumannBoundary = "ALL"
#########
if algo == "RI":
approx = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb)
else:
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.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], duality = False))
/ approx.estimatorNormEngine.norm(
approx.getRHS(k0s[j], duality = False)))
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(0)),
1.05*np.max(norm)*np.ones(approx.mus.shape, 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(0)),
4.*np.max(resApp)*np.ones(approx.mus.shape, 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/tests/hfengines/fracture.py b/tests/hfengines/fracture.py
index a7624a2..576f130 100644
--- a/tests/hfengines/fracture.py
+++ b/tests/hfengines/fracture.py
@@ -1,59 +1,58 @@
# 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 pytest
import numpy as np
import ufl
import fenics as fen
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.solver.fenics import affine_warping
@pytest.mark.xfail(reason = "no graphical interface")
def test_fracture():
mu0 = [45. ** .5, .6]
solver2 = MFE(mu0 = mu0, H = 1., L = .75, delta = .05, n = 20,
verbosity = 0)
uh2 = solver2.solve(mu0)[0]
solver1 = MFEND(mu0 = mu0, H = 1., L = .75, delta = .05, n = 20,
verbosity = 0)
uh1 = solver1.solve(mu0[0])[0]
L = mu0[1]
y = fen.SpatialCoordinate(solver1.V.mesh())[1]
warp1, warpI1 = affine_warping(solver1.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver1.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)
solver1.plotmesh([warp, warpI], show = False, figsize = (7, 7))
solver1.plot(uh1, [warp, warpI], what = 'REAL', show = False)
from matplotlib import pyplot as plt
plt.close('all')
assert np.isclose(
- solver1.norm(solver1.residual(uh1, mu0[0]), dual = True)[0],
- solver2.norm(solver2.residual(uh2, mu0), dual = True)[0],
+ solver1.norm(solver1.residual(mu0[0], uh1), dual = True)[0],
+ solver2.norm(solver2.residual(mu0, uh2), dual = True)[0],
atol = 1e-5)
assert np.isclose(solver1.norm(uh1 - uh2), 0., atol = 1e-6)
-
diff --git a/tests/hfengines/helmholtz_elasticity.py b/tests/hfengines/helmholtz_elasticity.py
index bba96cc..1ced33e 100644
--- a/tests/hfengines/helmholtz_elasticity.py
+++ b/tests/hfengines/helmholtz_elasticity.py
@@ -1,58 +1,58 @@
# 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.vector_linear_problem import (
LinearElasticityHelmholtzProblemEngine,
LinearElasticityHelmholtzProblemEngineDamped)
from rod_3d import rod3Dsolver
def test_helmholtz_elastic_rod():
solverBase = rod3Dsolver()
solver = LinearElasticityHelmholtzProblemEngine()
solver.V = solverBase.V
solver.lambda_ = solverBase.lambda_
solver.mu_ = solverBase.mu_
solver.forcingTerm = solverBase.forcingTerm
solver.DirichletBoundary = solverBase.DirichletBoundary
solver.NeumannBoundary = solverBase.NeumannBoundary
mu = 10
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 0.17836028624665373, rtol = 1e-5)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
8.070977e-07, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
- solver.norm(solver.residual(uh, mu, duality = False)[0],
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
+ solver.norm(solver.residual(mu, uh, duality = False)[0],
dual = True, duality = False), rtol = 1e-1)
def test_helmholtz_elastic_rod_damped():
solverBase = rod3Dsolver()
solver = LinearElasticityHelmholtzProblemEngineDamped()
solver.V = solverBase.V
solver.lambda_ = solverBase.lambda_
solver.mu_ = solverBase.mu_
solver.forcingTerm = solverBase.forcingTerm
solver.DirichletBoundary = solverBase.DirichletBoundary
solver.NeumannBoundary = solverBase.NeumannBoundary
solver.eta = 10
mu = 10
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 0.17646530119044376, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, 10)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(10, uh)[0], dual = True),
6.7057338e-07, rtol = 1e-1)
diff --git a/tests/hfengines/helmholtz_external.py b/tests/hfengines/helmholtz_external.py
index 3e74d16..c126c9f 100644
--- a/tests/hfengines/helmholtz_external.py
+++ b/tests/hfengines/helmholtz_external.py
@@ -1,69 +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 .
#
import pytest
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)[0]
assert np.isclose(solver.norm(uh), 19.9362, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
4.25056407e-13, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
- solver.norm(solver.residual(uh, mu, duality = False)[0],
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
+ solver.norm(solver.residual(mu, uh, duality = False)[0],
dual = True, duality = False), 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)[0]
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.setthAs(solver1.thAs)
solver2.setbs(solver1.bs)
+ solver2.setthbs(solver1.thbs)
uh2 = solver2.solve(mu)[0]
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
@pytest.mark.xfail(reason = "no graphical interface")
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)[0]
solver.plotmesh(show = False, figsize = (7, 7))
assert np.isclose(solver.norm(uh), 62.113, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
9.62989935e-13, rtol = 1e-1)
from matplotlib import pyplot as plt
plt.close('all')
diff --git a/tests/hfengines/helmholtz_internal.py b/tests/hfengines/helmholtz_internal.py
index d6fde8c..7414310 100644
--- a/tests/hfengines/helmholtz_internal.py
+++ b/tests/hfengines/helmholtz_internal.py
@@ -1,97 +1,97 @@
# 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 os, shutil
import numpy as np
from rrompy.hfengines.linear_problem import (
HelmholtzSquareBubbleDomainProblemEngine, HelmholtzSquareBubbleProblemEngine,
HelmholtzSquareTransmissionProblemEngine)
def test_helmholtz_square_io():
solver = HelmholtzSquareBubbleProblemEngine(kappa = 4, theta = 1., n = 20,
verbosity = 0)
mu = 5
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 145.0115, rtol = 1e-3)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
1.19934e-11, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
- solver.norm(solver.residual(uh, mu, duality = False)[0],
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
+ solver.norm(solver.residual(mu, uh, duality = False)[0],
dual = True, duality = False), rtol = 1e-1)
if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache")
filesOut = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".pvd" and x[:9] == "outSquare")]
filesOutData = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".vtu" and x[:9] == "outSquare")]
for fileOut in filesOut:
os.remove("./.pytest_cache/" + fileOut)
for fileOut in filesOutData:
os.remove("./.pytest_cache/" + fileOut)
solver.outParaview(uh, what = ["MESH", "ABS"],
filename = ".pytest_cache/outSquare",
forceNewFile = False)
filesOut = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".pvd" and x[:9] == "outSquare")]
filesOutData = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".vtu" and x[:9] == "outSquare")]
assert len(filesOut) == 1
assert len(filesOutData) == 1
os.remove("./.pytest_cache/" + filesOut[0])
os.remove("./.pytest_cache/" + filesOutData[0])
def test_helmholtz_transmission_io():
solver = HelmholtzSquareTransmissionProblemEngine(nT = 1, nB = 2,
theta = np.pi * 40 / 180, kappa = 4., n = 20, verbosity = 0)
mu = 5.
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 138.6609, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
3.7288565e-12, rtol = 1e-1)
if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache")
solver.outParaviewTimeDomain(uh, omega = mu,
filename = ".pytest_cache/outTrans",
forceNewFile = False, folder = True)
filesOut = [x for x in os.listdir("./.pytest_cache/outTrans") if
(x[-4:] == ".pvd" and x[:8] == "outTrans")]
filesOutData = [x for x in os.listdir("./.pytest_cache/outTrans") if
(x[-4:] == ".vtu" and x[:8] == "outTrans")]
assert len(filesOut) == 1
assert len(filesOutData) == 20
shutil.rmtree("./.pytest_cache/outTrans")
def test_helmholtz_domain_io():
solver = HelmholtzSquareBubbleDomainProblemEngine(kappa = 4, theta = 1.,
n = 10, mu0 = 1.5, verbosity = 0)
mu = 1.5
uh = solver.solve(mu)[0]
if not os.path.isdir("./.pytest_cache"): os.mkdir("./.pytest_cache")
solver.plot(uh, save = "./.pytest_cache/outDomain", show = False)
filesOut = [x for x in os.listdir("./.pytest_cache") if
(x[-4:] == ".eps" and x[:9] == "outDomain")]
assert len(filesOut) == 1
os.remove("./.pytest_cache/" + filesOut[0])
assert np.isclose(solver.norm(uh), 10.07843, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
6.14454989e-13, rtol = 1e-1)
diff --git a/tests/hfengines/laplace.py b/tests/hfengines/laplace.py
index 409ea74..0a3bed5 100644
--- a/tests/hfengines/laplace.py
+++ b/tests/hfengines/laplace.py
@@ -1,42 +1,42 @@
# 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 LaplaceDiskGaussian
from rrompy.hfengines.linear_problem.bidimensional import LaplaceDiskGaussian2
def test_laplace_disk():
solver = LaplaceDiskGaussian(n = 20, verbosity = 0)
mu = 1.5
solver.setSolver("BICG", {"tol" : 1e-15})
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 1.053403077447029, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
4.66728e-10, atol = 1e-7)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
- solver.norm(solver.residual(uh, mu, duality = False)[0],
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
+ solver.norm(solver.residual(mu, uh, duality = False)[0],
dual = True, duality = False), rtol = 1e-1)
def test_laplace_disk_2():
solver = LaplaceDiskGaussian2(n = 20, verbosity = 0)
mu = [[0., 1.5]]
mu = [0., 1.5]
uh = solver.solve(mu)[0]
assert np.isclose(solver.norm(uh), 1.0534030774205372, rtol = 1e-2)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
2.40043363e-08, atol = 1e-7)
diff --git a/tests/hfengines/linear_elasticity.py b/tests/hfengines/linear_elasticity.py
index 309716a..c4b294d 100644
--- a/tests/hfengines/linear_elasticity.py
+++ b/tests/hfengines/linear_elasticity.py
@@ -1,41 +1,41 @@
# 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.vector_linear_problem import (
LinearElasticityBeamPoissonRatio)
from rod_3d import rod3Dsolver
def test_elastic_beam():
solver = LinearElasticityBeamPoissonRatio(n = 10, rho_ = 1e3, g = 3,
E = 1e6, nu0 = .45, length = 5, verbosity = 0)
mu = .45
uh = solver.solve(mu)[0]
- assert np.isclose(solver.norm(uh), 5.866888537228743e-08, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(uh), 9.33957e-08, rtol = 1e-1)
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
8.4545952e-13, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
- solver.norm(solver.residual(uh, mu, duality = False)[0],
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
+ solver.norm(solver.residual(mu, uh, duality = False)[0],
dual = True, duality = False), rtol = 1e-1)
def test_elastic_rod():
solver = rod3Dsolver()
uh = solver.solve()[0]
assert np.isclose(solver.norm(uh), 0.15563476339534466, rtol = 1e-5)
- assert np.isclose(solver.norm(solver.residual(uh)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual([], uh)[0], dual = True),
1.2210129e-07, rtol = 1e-1)
diff --git a/tests/hfengines/matrix.py b/tests/hfengines/matrix.py
index d0959e5..d65dbd1 100644
--- a/tests/hfengines/matrix.py
+++ b/tests/hfengines/matrix.py
@@ -1,65 +1,57 @@
# 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 sp
-from rrompy.hfengines.base import MatrixEngineBase as MEB
+from rrompy.hfengines.base import MatrixEngineBase
def test_deterministic():
+ solver = MatrixEngineBase(verbosity = 0)
N = 100
- verb = 0
-
- solver = MEB(verbosity = verb)
solver.npar = 1
- solver.nAs = 2
- mu = 10. + .5j
-
+ solver.nAs, solver.nbs = 2, 1
solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- - sp.eye(N)]
- solver.nbs = 1
+ - sp.eye(N)]
solver.bs = [np.exp(1.j * np.linspace(0, -np.pi, N))]
+ solver.thAs = solver.getMonomialWeights(solver.nAs)
+ solver.thbs = solver.getMonomialWeights(solver.nbs)
+ mu = 10. + .5j
uh = solver.solve(mu)[0]
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
1.088e-15, rtol = 1e-1)
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
- solver.norm(solver.residual(uh, mu, duality = False)[0],
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
+ solver.norm(solver.residual(mu, uh, duality = False)[0],
dual = True, duality = False), rtol = 1e-1)
def test_random():
+ solver = MatrixEngineBase(verbosity = 0)
N = 100
- verb = 0
-
- solver = MEB(verbosity = verb)
+ solver.setSolver("SOLVE")
solver.npar = 1
- solver.nAs = 2
- mu = 1. + .5j
-
+ solver.nAs, solver.nbs = 2, 1
np.random.seed(420)
-
- 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)]
-
- solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
-
+ solver.thAs = solver.getMonomialWeights(solver.nAs)
+ solver.thbs = solver.getMonomialWeights(solver.nbs)
+ mu = 1. + .5j
uh = solver.solve(mu)[0]
-
- assert np.isclose(solver.norm(solver.residual(uh, mu)[0], dual = True),
+ assert np.isclose(solver.norm(solver.residual(mu, uh)[0], dual = True),
7.18658e-14, rtol = 1e-1)
diff --git a/tests/reduction_methods_1D/matrix_fft.py b/tests/reduction_methods_1D/matrix_fft.py
index d5ea4c0..d8e87b7 100644
--- a/tests/reduction_methods_1D/matrix_fft.py
+++ b/tests/reduction_methods_1D/matrix_fft.py
@@ -1,35 +1,37 @@
# 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.base import MatrixEngineBase as MEB
+from rrompy.hfengines.base import MatrixEngineBase
+
+class matrixFFT(MatrixEngineBase):
+ def __init__(self):
+ super().__init__(verbosity = 0)
+ N = 100
+ np.random.seed(420)
+ self.setSolver("SOLVE")
+ self.npar = 1
+ self.nAs, self.nbs = 2, 1
+ self.mu0 = 0.
+ fftB = np.fft.fft(np.eye(N)) * N**-.5
+ self.As = [fftB.dot(np.multiply(np.arange(1, 1 + N), fftB.conj()).T),
+ - np.eye(N)]
+ self.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
+ self.thAs = self.getMonomialWeights(self.nAs)
+ self.thbs = self.getMonomialWeights(self.nbs)
-def matrixFFT():
- N = 100
- solver = MEB(verbosity = 0)
- np.random.seed(420)
- solver.setSolver("SOLVE")
- fftB = np.fft.fft(np.eye(N)) * N**-.5
- solver.npar = 1
- solver.mu0 = 0.
- solver.nAs = 2
- solver.As = [fftB.dot(np.multiply(np.arange(1, 1 + N), fftB.conj()).T),
- - np.eye(N)]
- solver.nbs = 1
- solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
- return solver
diff --git a/tests/reduction_methods_1D/reduced_basis_1d.py b/tests/reduction_methods_1D/reduced_basis_1d.py
index 3cd2a78..715db16 100644
--- a/tests/reduction_methods_1D/reduced_basis_1d.py
+++ b/tests/reduction_methods_1D/reduced_basis_1d.py
@@ -1,57 +1,57 @@
# 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.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
from rrompy.parameter import checkParameterList
def test_LS():
solver = matrixFFT()
params = {"POD": True, "R": 5, "S": 10,
"sampler": QS([1., 7.], "CHEBYSHEV")}
approx = RB(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7)
approx.POD = False
approx.setupApprox()
for mu in approx.mus[approx.R :]:
assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-3)
-
+
def test_interp():
solver = matrixFFT()
params = {"POD": False, "S": 10, "sampler": QS([1., 7.], "CHEBYSHEV")}
approx = RB(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7)
def test_hermite():
mu = 1.5
solver = matrixFFT()
sampler0 = QS([1., 7.], "CHEBYSHEV")
points, _ = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
params = {"POD": True, "S": 12, "sampler": MS([1., 7.], points = points)}
approx = RB(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
diff --git a/tests/reduction_methods_1D/reduced_basis_greedy_1d.py b/tests/reduction_methods_1D/reduced_basis_greedy_1d.py
index 0086ab0..e56f2bb 100644
--- a/tests/reduction_methods_1D/reduced_basis_greedy_1d.py
+++ b/tests/reduction_methods_1D/reduced_basis_greedy_1d.py
@@ -1,59 +1,58 @@
# 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.greedy import ReducedBasisGreedy as RBG
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
def test_lax_tolerance(capsys):
mu = 2.25
solver = matrixFFT()
params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"),
"S": 4, "greedyTol": 1e-2,
"trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")}
approx = RBG(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 len(approx.mus) == 10
_, _, maxEst = approx.getMaxErrorEstimator(approx.muTest)
assert maxEst < 1e-2
assert np.isclose(approx.normErr(mu)[0], 1.5056e-05, rtol = 1e-1)
def test_samples_at_poles():
solver = matrixFFT()
params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"),
"S": 4, "nTestPoints": 100, "greedyTol": 1e-5,
"trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")}
approx = RBG(solver, 4., params, verbosity = 0)
approx.greedy()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0] / (1e-15+approx.normHF(mu)[0]),
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(np.array(approx.mus(0)) - lambda_)),
0., atol = 1e-1)
-
diff --git a/tests/reduction_methods_multiD/matrix_random.py b/tests/reduction_methods_multiD/matrix_random.py
index c82fad7..6765147 100644
--- a/tests/reduction_methods_multiD/matrix_random.py
+++ b/tests/reduction_methods_multiD/matrix_random.py
@@ -1,38 +1,40 @@
# 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.base import MatrixEngineBase as MEB
+from rrompy.hfengines.base import MatrixEngineBase
+
+class matrixRandom(MatrixEngineBase):
+ def __init__(self):
+ super().__init__(verbosity = 0)
+ N = 100
+ np.random.seed(420)
+ self.setSolver("SOLVE")
+ self.npar = 2
+ self.nAs, self.nbs = 3, 1
+ self.mu0 = [0., 0.]
+ d1 = np.random.randn(N)
+ Q1, _ = np.linalg.qr(np.random.randn(N, N))
+ d2 = np.random.randn(N)
+ Q2, _ = np.linalg.qr(np.random.randn(N, N))
+ self.As = [np.eye(N), Q1.dot(np.multiply(d1, Q1.conj()).T),
+ Q2.dot(np.multiply(d2, Q2.conj()).T)]
+ self.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
+ self.thAs = self.getMonomialWeights(self.nAs)
+ self.thbs = self.getMonomialWeights(self.nbs)
-def matrixRandom():
- N = 100
- solver = MEB(verbosity = 0)
- np.random.seed(420)
- solver.setSolver("SOLVE")
- solver.npar = 2
- solver.mu0 = [0., 0.]
- solver.nAs = 3
- d1 = np.random.randn(N)
- Q1, _ = np.linalg.qr(np.random.randn(N, N))
- d2 = np.random.randn(N)
- Q2, _ = np.linalg.qr(np.random.randn(N, N))
- solver.As = [np.eye(N), Q1.dot(np.multiply(d1, Q1.conj()).T),
- Q2.dot(np.multiply(d2, Q2.conj()).T)]
- solver.nbs = 1
- solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
- return solver
diff --git a/tests/utilities/poly_fitting.py b/tests/utilities/poly_fitting.py
index 936b051..909452e 100644
--- a/tests/utilities/poly_fitting.py
+++ b/tests/utilities/poly_fitting.py
@@ -1,116 +1,116 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
polydomcoeff, polyval,
- polyroots, polyvander,
- nextDerivativeIndices)
+ polyroots, polyvander)
+from rrompy.utilities.numerical import nextDerivativeIndices
from rrompy.parameter import checkParameterList
def test_chebyshev():
assert "CHEBYSHEV" in polybases
fitname = polyfitname("CHEBYSHEV")
domcoeff = polydomcoeff(5, "CHEBYSHEV")
assert fitname == "chebfit"
assert np.isclose(domcoeff, 16, rtol = 1e-5)
assert np.allclose(polyroots((-1, 1, -1, 1), "CHEBYSHEV"),
np.array([-.5, 0., 1.]), rtol = 1e-5)
Phi = polyvander(np.arange(5), 4, "CHEBYSHEV")
y = 2. * np.arange(5)
cFit = customFit(Phi, y)
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
assert np.allclose(polyval(np.arange(5), cFit, "CHEBYSHEV"), y,
rtol = 1e-5)
assert np.allclose(polyval(np.arange(5), cFit, "CHEBYSHEV", m = 1),
2. * np.ones(5), rtol = 1e-5)
def test_legendre():
assert "LEGENDRE" in polybases
fitname = polyfitname("LEGENDRE")
domcoeff = polydomcoeff([0, 5], "LEGENDRE")
assert fitname == "legfit"
assert np.allclose(domcoeff, [1., 63. / 8], rtol = 1e-5)
assert np.allclose(polyroots((1, 2, 3, 4), "LEGENDRE"),
np.array([-0.85099543, -0.11407192, 0.51506735]),
rtol = 1e-5)
Phi = polyvander(np.arange(5), 4, "LEGENDRE")
y = 2. * np.arange(5)
cFit = customFit(Phi, y)
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
assert np.allclose(polyval(np.arange(5), cFit, "LEGENDRE"), y, rtol = 1e-5)
assert np.allclose(polyval(np.arange(5), cFit, "LEGENDRE", m = 1),
2. * np.ones(5), rtol = 1e-5)
def test_monomial():
assert "MONOMIAL" in polybases
fitname = polyfitname("MONOMIAL")
domcoeff = polydomcoeff(5, "MONOMIAL")
assert fitname == "polyfit"
assert np.isclose(domcoeff, 1., rtol = 1e-5)
assert np.allclose(polyroots([0.+0.j, 1.+0.j, 0.+0.j, 1.+0.j], "MONOMIAL"),
np.array([0., 1.j, -1.j]), rtol = 1e-5)
Phi = polyvander(np.arange(5), 4, "MONOMIAL")
y = 2. * np.arange(5)
cFit = customFit(Phi, y)
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
assert np.allclose(polyval(np.arange(5), cFit, "MONOMIAL"), y, rtol = 1e-5)
assert np.allclose(polyval(np.arange(5), cFit, "MONOMIAL", m = 1),
2. * np.ones(5), rtol = 1e-5)
def test_cheb_confluence():
x = np.arange(5)
x = np.delete(x, 3)
Phi = polyvander(x, 4, "CHEBYSHEV", [[0, 1]] + [[0]] * 3,
reorder = [0, 2, 3, 1, 4])
y = 2. * np.arange(5)
y[3] = 2.
cFit = customFit(Phi, y)
mask = np.arange(len(y)) != 3
assert np.allclose(cFit, [0, 2, 0, 0, 0], rtol = 1e-5)
assert np.allclose(polyval(x, cFit, "CHEBYSHEV"), y[mask],
rtol = 1e-5)
assert np.allclose(polyval(x[0], cFit, "CHEBYSHEV", m = 1), y[~mask],
rtol = 1e-5)
def test_mon_confluence_2d():
x, _ = checkParameterList([[0, 0], [1, 1], [-1, -1]])
y = np.array([3., 5., 1., 2., 12., -2.]).reshape((-1, 1)) # 3+y+5x+2xy+x2y
Phi = polyvander(x, [2, 1], "MONOMIAL",
[[[0, 0], [1, 0], [0, 1], [1, 1]]] + [[[0, 0]]] * 2)
cFit = customFit(Phi, y).reshape((3, 2))
mask = np.array([0, 4, 5])
assert np.allclose(cFit.flatten(), [3, 1, 5, 2, 0, 1], atol = 1e-5)
assert np.allclose(polyval(x, cFit, "MONOMIAL").flatten(),
y[mask].flatten(), rtol = 1e-5)
assert np.allclose(polyval([x[0]], cFit, "MONOMIAL", m = [1, 0]), y[1],
rtol = 1e-5)
assert np.allclose(polyval([x[0]], cFit, "MONOMIAL", m = [0, 1]), y[2],
rtol = 1e-5)
assert np.allclose(polyval([x[0]], cFit, "MONOMIAL", m = [1, 1]), y[3],
rtol = 1e-5)
def test_derivative_indices_4d():
idxs = nextDerivativeIndices([], 4, 70)
idxMag = [np.sum(idx) for idx in idxs]
idxMagUnique, idxMagCount = np.unique(idxMag, return_counts = True)
idxMagCount = np.cumsum(idxMagCount)
assert np.allclose(idxMagUnique, np.arange(5), atol = 1e-10)
assert np.allclose(idxMagCount, [1, 5, 15, 35, 70], atol = 1e-10)
diff --git a/tests/utilities/radial_fitting.py b/tests/utilities/radial_fitting.py
index 7381c16..dd996f0 100644
--- a/tests/utilities/radial_fitting.py
+++ b/tests/utilities/radial_fitting.py
@@ -1,167 +1,167 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.poly_fitting.radial_basis import (radialGaussian,
thinPlateSpline,
multiQuadric,
polybases, polyfitname,
polydomcoeff,
polyval, polyvander,
polyvanderTotal)
-from rrompy.utilities.poly_fitting.polynomial import degreeTotalToFull
+from rrompy.utilities.numerical import degreeTotalToFull
from rrompy.parameter import checkParameterList
def test_monomial_gaussian():
polyrbname = "MONOMIAL_GAUSSIAN"
assert polyrbname in polybases
fitname = polyfitname(polyrbname)
domcoeff = polydomcoeff(5, polyrbname)
assert fitname == "polyfit_gaussian"
assert np.isclose(domcoeff, 1., rtol = 1e-5)
directionalWeights = np.array([5.])
xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
globalCoeffs = cRBCoeffs[4 :]
localCoeffs = cRBCoeffs[: 4]
ySupp = 1 + 2. * xSupp.data - .5 * xSupp.data ** 2.
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1)[0], globalCoeffs, localCoeffs,
xSupp, directionalWeights, polyrbname)
yyman = 1. + 2. * xx - .5 * xx ** 2.
for j, xc in enumerate(np.arange(-1, 3)):
r2j = (5. * (xx - xc)) ** 2.
rbj = radialGaussian(r2j)
assert np.allclose(rbj, np.exp(-.5 * r2j))
yyman += localCoeffs[j] * rbj
ySupp += localCoeffs[j] * radialGaussian((directionalWeights[0]
* (xSupp.data - xc)) ** 2.)
assert np.allclose(yy, yyman, atol = 1e-5)
VanT = polyvander(xSupp, [2], polyrbname,
directionalWeights = directionalWeights)
ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
out = customFit(VanT, ySupp)
assert np.allclose(out, cRBCoeffs, atol = 1e-8)
def test_legendre_thinplate():
polyrbname = "LEGENDRE_THINPLATE"
assert polyrbname in polybases
fitname = polyfitname(polyrbname)
domcoeff = polydomcoeff(5, polyrbname)
assert fitname == "legfit_thinplate"
assert np.isclose(domcoeff, 63. / 8, rtol = 1e-5)
directionalWeights = np.array([.5])
xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
localCoeffs = cRBCoeffs[: 4]
globalCoeffs = cRBCoeffs[4 :]
ySupp = 1 + 2. * xSupp.data - .5 * (.5 * (3. * xSupp.data ** 2. - 1.))
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1)[0], globalCoeffs, localCoeffs,
xSupp, directionalWeights, polyrbname)
yyman = 1. + 2. * xx - .5 * (.5 * (3. * xx ** 2. - 1.))
for j, xc in enumerate(np.arange(-1, 3)):
r2j = (directionalWeights[0] * (xx - xc)) ** 2.
rbj = thinPlateSpline(r2j)
assert np.allclose(rbj, .5 * r2j * np.log(np.finfo(float).eps + r2j))
yyman += localCoeffs[j] * rbj
ySupp += localCoeffs[j] * thinPlateSpline((directionalWeights[0]
* (xSupp.data - xc)) ** 2.)
assert np.allclose(yy, yyman, atol = 1e-5)
VanT = polyvander(xSupp, [2], polyrbname,
directionalWeights = directionalWeights)
ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
out = customFit(VanT, ySupp)
assert np.allclose(out, cRBCoeffs, atol = 1e-8)
def test_chebyshev_multiquadric():
polyrbname = "CHEBYSHEV_MULTIQUADRIC"
assert polyrbname in polybases
fitname = polyfitname(polyrbname)
domcoeff = polydomcoeff(5, polyrbname)
assert fitname == "chebfit_multiquadric"
assert np.isclose(domcoeff, 16, rtol = 1e-5)
directionalWeights = np.array([1.])
xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
localCoeffs = cRBCoeffs[: 4]
globalCoeffs = cRBCoeffs[4 :]
ySupp = 1 + 2. * xSupp.data - .5 * (2. * xSupp.data ** 2. - 1.)
xx = np.linspace(-2., 3., 100)
yy = polyval(checkParameterList(xx, 1)[0], globalCoeffs, localCoeffs,
xSupp, directionalWeights, polyrbname)
yyman = 1. + 2. * xx - .5 * (2. * xx ** 2. - 1.)
for j, xc in enumerate(np.arange(-1, 3)):
r2j = (directionalWeights[0] * (xx - xc)) ** 2.
rbj = multiQuadric(r2j)
assert np.allclose(rbj, np.power(r2j + 1, -.5))
yyman += localCoeffs[j] * rbj
ySupp += localCoeffs[j] * multiQuadric((directionalWeights[0]
* (xSupp.data - xc)) ** 2.)
assert np.allclose(yy, yyman, atol = 1e-5)
VanT = polyvander(xSupp, [2], polyrbname,
directionalWeights = directionalWeights)
ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
out = customFit(VanT, ySupp)
assert np.allclose(out, cRBCoeffs, atol = 1e-8)
def test_total_degree_2d():
values = lambda x, y: (x - 3.) ** 2. * y - (x + 1.) * y ** 2.
polyrbname = "CHEBYSHEV_GAUSSIAN"
xs, ys = np.meshgrid(np.linspace(0., 4., 5), np.linspace(0., 4., 4))
xySupp = np.concatenate((xs.reshape(-1, 1), ys.reshape(-1, 1)), axis = 1)
zs = values(xs, ys)
zSupp = zs.flatten()
deg = 3
directionalWeights = [2., 1.]
VanT, _, reidxs = polyvanderTotal(xySupp, deg, polyrbname,
directionalWeights = directionalWeights)
VanT = VanT[reidxs]
VanT = VanT[:, reidxs]
cFit = np.linalg.solve(VanT, np.pad(zSupp, (0, len(VanT) - len(zSupp)),
'constant'))
globCoeff = degreeTotalToFull([deg + 1] * 2, 2, cFit[len(zSupp) :])
localCoeffs = cFit[: len(zSupp)]
globalCoeffs = globCoeff
xx, yy = np.meshgrid(np.linspace(0., 4., 100), np.linspace(0., 4., 100))
xxyy = np.concatenate((xx.reshape(-1, 1), yy.reshape(-1, 1)), axis = 1)
zz = polyval(xxyy, globalCoeffs, localCoeffs, xySupp, directionalWeights,
polyrbname).reshape(xx.shape)
zzex = values(xx, yy)
error = np.abs(zz - zzex)
print(np.max(error))
assert np.max(error) < 1e-10
diff --git a/tests/utilities/sampling.py b/tests/utilities/sampling.py
index a441267..4df35a7 100644
--- a/tests/utilities/sampling.py
+++ b/tests/utilities/sampling.py
@@ -1,78 +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
import scipy.sparse as sp
-from rrompy.hfengines.base import MatrixEngineBase as MEB
+from rrompy.hfengines.base import MatrixEngineBase
from rrompy.sampling.standard import (SamplingEngineStandard,
SamplingEngineStandardPOD)
from rrompy.parameter import parameterList
+class matrixEngine(MatrixEngineBase):
+ def __init__(self):
+ super().__init__(verbosity = 0)
+ N = 100
+ self.npar = 1
+ self.nAs, self.nbs = 2, 1
+ self.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
+ - sp.eye(N)]
+ self.bs = [np.exp(1.j * np.linspace(0, -np.pi, N))]
+ self.thAs = self.getMonomialWeights(self.nAs)
+ self.thbs = self.getMonomialWeights(self.nbs)
+
def test_krylov():
- N = 100
mu = 10. + .5j
- solver = MEB(verbosity = 0)
- solver.npar = 1
- solver.nAs = 2
-
- solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- - sp.eye(N)]
- solver.nbs = 1
- solver.bs = [np.exp(1.j * np.linspace(0, -np.pi, N))]
+ solver = matrixEngine()
samplingEngine = SamplingEngineStandard(solver, verbosity = 0)
-
samples = samplingEngine.iterSample([mu] * 5).data
assert samples.shape == (100, 5)
assert np.isclose(np.linalg.norm(samples), 37.02294804524299, rtol = 1e-5)
def test_distributed():
- N = 100
mus = parameterList(np.linspace(5, 15, 11) + .5j)
- solver = MEB(verbosity = 0)
- solver.npar = 1
- solver.nAs = 2
-
- solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- - sp.eye(N)]
- solver.nbs = 1
- solver.bs = [np.exp(1.j * np.linspace(0, -np.pi, N))]
+ solver = matrixEngine()
samplingEngine = SamplingEngineStandard(solver, verbosity = 0)
-
samples = samplingEngine.iterSample(mus).data
assert samples.shape == (100, 11)
assert np.isclose(np.linalg.norm(samples), 8.59778606421386, rtol = 1e-5)
def test_distributed_pod():
- N = 100
mus = np.linspace(5, 15, 11) + .5j
- solver = MEB(verbosity = 0)
- solver.npar = 1
- solver.nAs = 2
-
- solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- - sp.eye(N)]
- solver.nbs = 1
- solver.bs = [np.exp(1.j * np.linspace(0, -np.pi, N))]
+ solver = matrixEngine()
samplingEngine = SamplingEngineStandardPOD(solver, verbosity = 0)
samples = samplingEngine.iterSample(mus).data
assert samples.shape == (100, 11)
assert np.isclose(np.linalg.norm(samples), 3.3166247903553994, rtol = 1e-5)
assert np.isclose(np.linalg.cond(samples.conj().T.dot(samples)), 1.,
rtol = 1e-5)