diff --git a/examples/2d/base/matrix_passive.py b/examples/2d/base/matrix_passive.py new file mode 100644 index 0000000..cc1331a --- /dev/null +++ b/examples/2d/base/matrix_passive.py @@ -0,0 +1,13 @@ +from rrompy.hfengines.linear_problem.tridimensional import \ + MatrixDynamicalPassive as MDP + +verb = 100 +n = 100 +b = 10 + +mu0 = [0., 10] +mutar = [4., 15] + +solver = MDP(mu0 = mu0, n = n, b = b, verbosity = verb) +uh = solver.solve(mutar)[0] +solver.plot(uh) diff --git a/examples/2d/pod/fracture_pod.py b/examples/2d/pod/fracture_pod.py index 3ab7db4..957eb08 100644 --- a/examples/2d/pod/fracture_pod.py +++ b/examples/2d/pod/fracture_pod.py @@ -1,218 +1,277 @@ import numpy as np from rrompy.hfengines.linear_problem.bidimensional import \ MembraneFractureEngine as MFE from rrompy.reduction_methods.distributed import RationalInterpolant as RI from rrompy.reduction_methods.distributed import RBDistributed as RBD +from rrompy.reduction_methods.pivoting import RBDistributedPivoted as RBP +from rrompy.reduction_methods.pivoting import \ + RationalInterpolantPivotedPoleMatching as RIPPM from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, QuadratureSamplerTotal as QST, ManualSampler as MS, RandomSampler as RS) verb = 5 -size = 1 +size = 2 show_sample = False -show_norm = False +show_norm = True clip = -1 #clip = .4 #clip = .6 -homogeneize = False -#homogeneize = True Delta = 0 MN = 5 R = (MN + 2) * (MN + 1) // 2 S = [int(np.ceil(R ** .5))] * 2 PODTol = 1e-8 +SPivot = [MN + 1, 4] +MMarginal = SPivot[1] - 2 +matchingWeight = 10. samples = "centered" samples = "distributed" +samples = "pivoted" algo = "rational" -#algo = "RB" +algo = "RB" sampling = "quadrature" sampling = "quadrature_total" #sampling = "random" +samplingM = "quadrature" +#samplingM = "quadrature_total" +#samplingM = "random" if samples == "distributed": radial = 0 # radial = "gaussian" # radial = "thinplate" # radial = "multiquadric" rW0 = 5. radialWeight = [rW0] * 2 +if samples == "pivoted": + radial = 0 +# radial = "gaussian" +# radial = "thinplate" +# radial = "multiquadric" + rW0 = 5. + radialWeight = [rW0] + radialM = 0 + 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) rescaling = [lambda x: np.power(x, 2.), lambda x: x] rescalingInv = [lambda x: np.power(x, .5), lambda x: x] if algo == "rational": params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True} if samples == "distributed": params['polybasis'] = "CHEBYSHEV" # params['polybasis'] = "LEGENDRE" # params['polybasis'] = "MONOMIAL" params['E'] = MN params['radialBasis'] = radial params['radialBasisWeights'] = radialWeight method = RI elif samples == "centered": params['polybasis'] = "MONOMIAL" params['S'] = R method = RI + elif samples == "pivoted": + params['S'] = [SPivot[0]] + params['SMarginal'] = [SPivot[1]] + params['MMarginal'] = MMarginal + params['polybasisPivot'] = "CHEBYSHEV" + params['polybasisMarginal'] = "MONOMIAL" + params['matchingWeight'] = matchingWeight + params['radialBasisPivot'] = radial + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsPivot'] = radialWeight + params['radialBasisWeightsMarginal'] = radialWeightM + method = RIPPM else: #if algo == "RB": params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S, 'POD':True, 'PODTolerance':PODTol} if samples == "distributed": method = RBD elif samples == "centered": params['S'] = R method = RBD + elif samples == "pivoted": + params['S'] = [SPivot[0]] + params['SMarginal'] = [SPivot[1]] + params['MMarginal'] = MMarginal + params['polybasisMarginal'] = "MONOMIAL" + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsMarginal'] = radialWeightM + method = RBP if samples == "distributed": if sampling == "quadrature": params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) # params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling, # scalingInv = rescalingInv) # params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling, # scalingInv = rescalingInv) params['S'] = [max(j, MN + 1) for j in params['S']] elif sampling == "quadrature_total": params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R else: # if sampling == "random": params['sampler'] = RS(murange, "HALTON", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R - elif samples == "centered": params['sampler'] = MS(murange, points = [mu0], scaling = rescaling, scalingInv = rescalingInv) +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") + -approx = method(solver, mu0 = mu0, approxParameters = params, - verbosity = verb, homogeneized = homogeneize) -if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False +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', homogeneized = False, what = "REAL") approx.plotHF(mutar, [warp, warpI], name = 'u_HF', homogeneized = False, what = "REAL") approx.plotErr(mutar, [warp, warpI], name = 'err', homogeneized = False, what = "REAL") # approx.plotRes(mutar, [warp, warpI], name = 'res', # homogeneized = False, what = "REAL") - appErr = approx.normErr(mutar, homogeneized = homogeneize) - solNorm = approx.normHF(mutar, homogeneized = homogeneize) - resNorm = approx.normRes(mutar, homogeneized = homogeneize) - RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize) + appErr = approx.normErr(mutar) + solNorm = approx.normHF(mutar) + resNorm = approx.normRes(mutar) + RHSNorm = approx.normRHS(mutar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) -if algo == "rational" and approx.N > 0: - from plot_zero_set import plotZeroSet2 - muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0, - 200, [2., 1.], clip = clip) +from plot_zero_set import plotZeroSet2 +muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0, + 200, [2., 1.], clip = clip) if show_norm: solver._solveBatchSize = 100 from plot_inf_set import plotInfSet2 muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2( - murange, murangeEff, approx, mu0, 50, - [2., 1.], clip = clip, relative = False, - normalizeDen = True) -print(approx.getPoles([None, .5]) ** 2.) + murange, murangeEff, approx, mu0, 25, + [2., 1.], clip = clip, relative = True, + nobeta = True) +print(np.sort(approx.getPoles([None, .5]) ** 2.)) diff --git a/examples/2d/pod/matrix_passive_pod.py b/examples/2d/pod/matrix_passive_pod.py new file mode 100644 index 0000000..9ccfd7b --- /dev/null +++ b/examples/2d/pod/matrix_passive_pod.py @@ -0,0 +1,197 @@ +import numpy as np +from rrompy.hfengines.linear_problem.tridimensional import \ + MatrixDynamicalPassive as MDP +from rrompy.reduction_methods.distributed import RationalInterpolant as RI +from rrompy.reduction_methods.distributed import RBDistributed as RBD +from rrompy.reduction_methods.pivoting import RationalInterpolantPivoted as RIP +from rrompy.reduction_methods.pivoting import RBDistributedPivoted as RBP +from rrompy.reduction_methods.pivoting import \ + RationalInterpolantPivotedPoleMatching as RIPPM +from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, + QuadratureSamplerTotal as QST, + ManualSampler as MS, + RandomSampler as RS) + +verb = 5 +size = 3 +show_sample = True +show_norm = True + +Delta = 0 +MN = 7 +R = (MN + 2) * (MN + 1) // 2 +S = [int(np.ceil(R ** .5))] * 2 +PODTol = 1e-6 +SPivot = [MN + 1, 3] +MMarginal = SPivot[1] - 1 + +samples = "centered" +samples = "distributed" +samples = "pivoted" +samples = "pole matching" +algo = "rational" +#algo = "RB" +sampling = "quadrature" +#sampling = "quadrature_total" +#sampling = "random" +samplingM = "quadrature" +#samplingM = "quadrature_total" +#samplingM = "random" + +if samples in ["distributed", "pivoted", "pole matching"]: + radial = 0 +# radial = "gaussian" +# radial = "thinplate" +# radial = "multiquadric" + rW0 = 10. + radialWeight = [rW0] +if samples in ["pivoted", "pole matching"]: + radialM = 0 + radialM = "gaussian" +# radialM = "thinplate" +# radialM = "multiquadric" + rMW0 = 2. + radialWeightM = [rMW0] +matchingWeight = 1. + +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':S, 'POD':True} + if samples == "distributed": + params['polybasis'] = "CHEBYSHEV" +# params['polybasis'] = "LEGENDRE" +# params['polybasis'] = "MONOMIAL" + params['E'] = MN + params['radialBasis'] = radial + params['radialBasisWeights'] = radialWeight + method = RI + elif samples == "centered": + params['polybasis'] = "MONOMIAL" + params['S'] = R + method = RI + elif samples in ["pivoted", "pole matching"]: + params['S'] = [SPivot[0]] + params['SMarginal'] = [SPivot[1]] + params['MMarginal'] = MMarginal + params['polybasisPivot'] = "CHEBYSHEV" + params['polybasisMarginal'] = "MONOMIAL" + params['radialBasisPivot'] = radial + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsPivot'] = radialWeight + params['radialBasisWeightsMarginal'] = radialWeightM + if samples == "pivoted": + method = RIP + else: + params['matchingWeight'] = matchingWeight + method = RIPPM +else: #if algo == "RB": + params = {'R':(MN + 2 + Delta) * (MN + 1 + Delta) // 2, 'S':S, + 'POD':True, 'PODTolerance':PODTol} + if samples == "distributed": + method = RBD + elif samples == "centered": + params['S'] = R + method = RBD + elif samples == "pivoted": + params['R'] = SPivot[0] + params['S'] = [SPivot[0]] + params['SMarginal'] = [SPivot[1]] + params['MMarginal'] = SPivot[1] - 2 + params['polybasisMarginal'] = "MONOMIAL" + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsMarginal'] = radialWeightM + method = RBP + +if samples == "distributed": + if sampling == "quadrature": + params['sampler'] = QS(murange, "CHEBYSHEV") +# params['sampler'] = QS(murange, "GAUSSLEGENDRE") +# params['sampler'] = QS(murange, "UNIFORM") + 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 in ["pivoted", "pole matching"]: + 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 not in ["pivoted", "pole matching"]: + 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', homogeneized = False, + what = "REAL") + approx.plotHF(mutar, name = 'u_HF', homogeneized = False, what = "REAL") + approx.plotErr(mutar, name = 'err', homogeneized = False, what = "REAL") +# approx.plotRes(mutar, name = 'res', homogeneized = False, what = "REAL") + appErr = approx.normErr(mutar) + solNorm = approx.normHF(mutar) + resNorm = approx.normRes(mutar) + RHSNorm = approx.normRHS(mutar) + print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, + np.divide(appErr, solNorm))) + print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, + np.divide(resNorm, RHSNorm))) + +if hasattr(approx.trainedModel, "getQVal") and approx.N > 0: + from plot_zero_set import plotZeroSet2 + muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0, + 200, [1., 1]) + +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 5f5c891..d7ba227 100644 --- a/examples/2d/pod/plot_inf_set.py +++ b/examples/2d/pod/plot_inf_set.py @@ -1,335 +1,335 @@ 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(mu2x, [ZTmin] * len(mu2x), 'kx') + plt.plot(np.real(mu2x), [ZTmin] * len(mu2x), 'kx') plt.plot([murangeExp[0][0]] * 2, [ZTmin, ZTmax], 'm:') plt.plot([murangeExp[1][0]] * 2, [ZTmin, ZTmax], 'm:') plt.xlim(mu1[0], mu1[-1]) plt.title("|u(mu)|, |u_app(mu)|") plt.grid() plt.show() plt.figure(figsize = (15, 7)) plt.jet() plt.semilogy(mu1, R) for l_ in approx.trainedModel.getPoles(): plt.plot([np.real(l_ ** exp)] * 2, [Rmin, Rmax], 'b:') - plt.plot(mu2x, [Rmax] * len(mu2x), 'kx') + plt.plot(np.real(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(): plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:') - plt.plot(mu2x, [Emax] * len(mu2x), 'kx') + plt.plot(np.real(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(): plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:') - plt.plot(mu2x, [betamax] * len(mu2x), 'kx') + plt.plot(np.real(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(): plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:') - plt.plot(mu2x, [Emax] * len(mu2x), 'kx') + plt.plot(np.real(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(): plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:') - plt.plot(mu2x, [etamax] * len(mu2x), 'kx') + plt.plot(np.real(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 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) 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) 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(np.real(mu2x), np.real(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(np.real(mu2x), np.real(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(np.real(mu2x), np.real(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(np.real(mu2x), np.real(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(np.real(mu2x), np.real(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(np.real(mu2x), np.real(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(np.real(mu2x), np.real(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 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)\ .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 69b32ed..825bf66 100644 --- a/examples/2d/pod/plot_zero_set.py +++ b/examples/2d/pod/plot_zero_set.py @@ -1,80 +1,93 @@ import warnings import numpy as np from matplotlib import pyplot as plt def plotZeroSet1(murange, murangeEff, approx, mu0, nSamples = 200, exp = 2.): if hasattr(approx, "mus"): mu2x = approx.mus(0) ** exp else: mu2x = mu0[0] ** exp murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]] mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp, nSamples) mus = np.power(mu1, 1. / exp) mu1 = np.real(mu1) Z = approx.trainedModel.getQVal(mus) Zabs = np.abs(Z) Zmin, Zmax = np.min(Zabs), np.max(Zabs) plt.figure(figsize = (15, 7)) plt.jet() plt.semilogy(mu1, Zabs) for l_ in approx.trainedModel.getPoles(): plt.plot([np.real(l_ ** exp)] * 2, [Zmin, Zmax], 'b--') plt.plot(mu2x, [Zmax] * len(mu2x), 'kx') plt.plot([murangeExp[0][0]] * 2, [Zmin, Zmax], 'm:') plt.plot([murangeExp[1][0]] * 2, [Zmin, Zmax], 'm:') plt.xlim(mu1[0], mu1[-1]) plt.title("|Q(mu)|") plt.grid() plt.show() return mus, Z def plotZeroSet2(murange, murangeEff, approx, mu0, nSamples = 200, - exps = [2., 2.], clip = -1): + 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] - 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: + 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) + try: + 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) plt.figure(figsize = (15, 7)) plt.jet() - p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap, - levels = np.linspace(Zabsmin, Zabsmax, 50)) - if clip > 0: - plt.contour(Mu1, Mu2, Zabs, [Zabsmin]) - plt.contour(Mu1, Mu2, np.real(Z), [0.], linestyles = 'dashed') - plt.contour(Mu1, Mu2, np.imag(Z), [0.], linewidths = 1, - linestyles = 'dotted') - plt.plot(mu2x, mu2y, 'kx') + 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') 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) + 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 5668df0..4cf18a2 100644 --- a/examples/2d/pod/square_pod.py +++ b/examples/2d/pod/square_pod.py @@ -1,143 +1,190 @@ import numpy as np from rrompy.hfengines.linear_problem.bidimensional import \ HelmholtzSquareDomainProblemEngine as HSDPE from rrompy.reduction_methods.distributed import RationalInterpolant as RI from rrompy.reduction_methods.distributed import RBDistributed as RBD +from rrompy.reduction_methods.pivoting import \ + RationalInterpolantPivotedPoleMatching as RIPPM from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, QuadratureSamplerTotal as QST, ManualSampler as MS, RandomSampler as RS) verb = 5 -size = 1 +size = 7 show_sample = False show_norm = True ignore_forcing = True ignore_forcing = False clip = -1 #clip = .4 #clip = .6 MN = 6 R = (MN + 2) * (MN + 1) // 2 S = [int(np.ceil(R ** .5))] * 2 +SPivot = [MN + 1, 3] +MMarginal = SPivot[1] - 2 +matchingWeight = 10. samples = "centered" samples = "distributed" +samples = "pivoted" algo = "rational" #algo = "RB" sampling = "quadrature" sampling = "quadrature_total" #sampling = "random" +samplingM = "quadrature" +#samplingM = "quadrature_total" +#samplingM = "random" + +if samples == "pivoted": + radialM = 0 + radialM = "gaussian" +# radialM = "thinplate" +# radialM = "multiquadric" + rMW0 = 2. + radialWeightM = [rMW0] if size == 1: # small mu0 = [4 ** .5, 1.5 ** .5] mutar = [5 ** .5, 1.75 ** .5] murange = [[2 ** .5, 1. ** .5], [6 ** .5, 2. ** .5]] elif size == 2: # medium mu0 = [4 ** .5, 1.75 ** .5] mutar = [5 ** .5, 1.25 ** .5] murange = [[1 ** .5, 1. ** .5], [7 ** .5, 2.5 ** .5]] elif size == 3: # fat mu0 = [6 ** .5, 4 ** .5] mutar = [2 ** .5, 2.5 ** .5] murange = [[0 ** .5, 2 ** .5], [12 ** .5, 6 ** .5]] elif size == 4: # crowded mu0 = [10 ** .5, 2 ** .5] mutar = [9 ** .5, 2.25 ** .5] murange = [[8 ** .5, 1.5 ** .5], [12 ** .5, 2.5 ** .5]] elif size == 5: # tall mu0 = [11 ** .5, 2.25 ** .5] mutar = [10.5 ** .5, 2.5 ** .5] murange = [[10 ** .5, 1.5 ** .5], [12 ** .5, 3 ** .5]] elif size == 6: # taller mu0 = [11 ** .5, 2.25 ** .5] mutar = [10.5 ** .5, 2.5 ** .5] murange = [[10 ** .5, 1.25 ** .5], [12 ** .5, 3.25 ** .5]] elif size == 7: # low mu0 = [7 ** .5, .75 ** .5] mutar = [6.5 ** .5, .9 ** .5] murange = [[6 ** .5, .5 ** .5], [8 ** .5, 1. ** .5]] -aEff = 1.25 +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 rescaling = [lambda x: np.power(x, 2.), lambda x: x] rescalingInv = [lambda x: np.power(x, .5), lambda x: x] if algo == "rational": params = {'N':MN, 'M':MN, 'S':S, 'POD':True} if samples == "distributed": params['polybasis'] = "CHEBYSHEV" # params['polybasis'] = "LEGENDRE" # params['polybasis'] = "MONOMIAL" params['E'] = MN method = RI elif samples == "centered": params['polybasis'] = "MONOMIAL" params['S'] = R method = RI + elif samples == "pivoted": + params['S'] = [SPivot[0]] + params['SMarginal'] = [SPivot[1]] + params['MMarginal'] = MMarginal + params['polybasisPivot'] = "CHEBYSHEV" + params['polybasisMarginal'] = "MONOMIAL" + params['matchingWeight'] = matchingWeight + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsMarginal'] = radialWeightM + method = RIPPM else: #if algo == "RB": params = {'R':R, 'S':S, 'POD':True} if samples == "distributed": method = RBD elif samples == "centered": params['S'] = R method = RBD if samples == "distributed": if sampling == "quadrature": params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) # params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling, # scalingInv = rescalingInv) # params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling, # scalingInv = rescalingInv) params['S'] = [max(j, MN + 1) for j in params['S']] elif sampling == "quadrature_total": params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R else: # if sampling == "random": params['sampler'] = RS(murange, "HALTON", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R - elif samples == "centered": params['sampler'] = MS(murange, points = [mu0], scaling = rescaling, scalingInv = rescalingInv) +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") -approx = method(solver, mu0 = mu0, approxParameters = params, - verbosity = verb) +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, 25, [2., 1.], clip = clip) + murange, murangeEff, approx, mu0, 20, [2., 1.], clip = clip, + relative = False, nobeta = True) diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py index b1b9a60..cba9424 100644 --- a/examples/3d/pod/fracture3_pod.py +++ b/examples/3d/pod/fracture3_pod.py @@ -1,195 +1,260 @@ 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.distributed import RationalInterpolant as RI from rrompy.reduction_methods.distributed import RBDistributed as RBD +from rrompy.reduction_methods.pivoting import \ + RationalInterpolantPivotedPoleMatching as RIPPM +from rrompy.reduction_methods.pivoting import \ + RBDistributedPivotedPoleMatching as RBP from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, QuadratureSamplerTotal as QST, ManualSampler as MS, RandomSampler as RS) -verb = 50 +verb = 70 size = 4 show_sample = False show_norm = False clip = -1 #clip = .4 #clip = .6 -homogeneize = False -#homogeneize = True Delta = 0 MN = 5 R = (MN + 3) * (MN + 2) * (MN + 1) // 6 S = [int(np.ceil(R ** (1. / 3.)))] * 3 PODTol = 1e-8 +SPivot = [MN + 1, 2, 2] +MMarginal = SPivot[1] - 1 +matchingWeight = 10. samples = "centered" samples = "distributed" +samples = "pivoted" algo = "rational" #algo = "RB" sampling = "quadrature" sampling = "quadrature_total" #sampling = "random" +samplingM = "quadrature" +#samplingM = "quadrature_total" +#samplingM = "random" if samples == "distributed": radial = 0 radial = "gaussian" # radial = "thinplate" # radial = "multiquadric" rW0 = 5. radialWeight = [rW0] * 3 +if samples == "pivoted": + radial = 0 +# radial = "gaussian" +# radial = "thinplate" +# radial = "multiquadric" + rW0 = 5. + radialWeight = [rW0] + radialM = 0 + 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]] 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) rescaling = [lambda x: np.power(x, 2.), lambda x: x, lambda x: x] rescalingInv = [lambda x: np.power(x, .5), lambda x: x, lambda x: x] if algo == "rational": params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True} if samples == "distributed": params['polybasis'] = "CHEBYSHEV" # params['polybasis'] = "LEGENDRE" # params['polybasis'] = "MONOMIAL" params['E'] = MN params['radialBasis'] = radial params['radialBasisWeights'] = radialWeight method = RI elif samples == "centered": params['polybasis'] = "MONOMIAL" params['S'] = R method = RI + elif samples == "pivoted": + params['S'] = [SPivot[0]] + params['SMarginal'] = SPivot[1:] + params['MMarginal'] = MMarginal + params['polybasisPivot'] = "CHEBYSHEV" + params['polybasisMarginal'] = "MONOMIAL" + params['matchingWeight'] = matchingWeight + params['radialBasisPivot'] = radial + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsPivot'] = radialWeight + params['radialBasisWeightsMarginal'] = radialWeightM + method = RIPPM else: #if algo == "RB": params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6, 'S':S, 'POD':True, 'PODTolerance':PODTol} if samples == "distributed": method = RBD elif samples == "centered": params['S'] = R method = RBD + elif samples == "pivoted": + params['S'] = [SPivot[0]] + params['SMarginal'] = SPivot[1:] + params['MMarginal'] = MMarginal + params['polybasisMarginal'] = "MONOMIAL" + params['matchingWeight'] = matchingWeight + params['radialBasisMarginal'] = radialM + params['radialBasisWeightsMarginal'] = radialWeightM + method = RBP if samples == "distributed": if sampling == "quadrature": params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) # params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling, # scalingInv = rescalingInv) # params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling, # scalingInv = rescalingInv) params['S'] = [max(j, MN + 1) for j in params['S']] elif sampling == "quadrature_total": params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R else: # if sampling == "random": params['sampler'] = RS(murange, "HALTON", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R elif samples == "centered": params['sampler'] = MS(murange, points = [mu0], scaling = rescaling, scalingInv = rescalingInv) +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") + params['S'] = MN + 1 + else: # if sampling == "random": + params['samplerPivot'] = RS([murange[0][0], murange[1][0]], "HALTON") + params['S'] = MN + 1 + 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") + params['SMarginal'] = (MMarginal + 2) * (MMarginal + 1) // 2 + params['polybasisMarginal'] = "CHEBYSHEV" + else: # if samplingM == "random": + params['samplerMarginal'] = RS([murange[0][1:], murange[1][1:]], "HALTON") + params['SMarginal'] = (MMarginal + 2) * (MMarginal + 1) // 2 -approx = method(solver, mu0 = mu0, approxParameters = params, - verbosity = verb, homogeneized = homogeneize) -if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False +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', homogeneized = False, what = "REAL") approx.plotHF(mutar, warps, name = 'u_HF', homogeneized = False, what = "REAL") approx.plotErr(mutar, warps, name = 'err', homogeneized = False, what = "REAL") # approx.plotRes(mutar, warps, name = 'res', # homogeneized = False, what = "REAL") - appErr = approx.normErr(mutar, homogeneized = homogeneize) - solNorm = approx.normHF(mutar, homogeneized = homogeneize) - resNorm = approx.normRes(mutar, homogeneized = homogeneize) - RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize) + 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(approx.trainedModel.data.mus(0) ** 2., - approx.trainedModel.data.mus(1), - approx.trainedModel.data.mus(2), '.') +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 -if algo == "rational" and approx.N > 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) +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_zero_set_3.py b/examples/3d/pod/plot_zero_set_3.py index 19b2856..c1db350 100644 --- a/examples/3d/pod/plot_zero_set_3.py +++ b/examples/3d/pod/plot_zero_set_3.py @@ -1,167 +1,170 @@ import warnings from copy import deepcopy as copy import numpy as np from matplotlib import pyplot as plt def plotZeroSet3(murange, murangeEff, approx, mu0, nSamples = 200, marginal = [None, None, .05], exps = [2., 1., 1.], - clip = -1, imagTol = 1e-8): + clip = -1, polesImTol : float = 1e-5): mNone = [i for i, m in enumerate(marginal) if m is None] nNone = len(mNone) if nNone not in [1, 2, 3]: raise if nNone == 1: exp = exps[mNone[0]] 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.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)] mu1 = np.real(mu1) Z = approx.trainedModel.getQVal(mus) Zabs = np.abs(Z) Zmin, Zmax = np.min(Zabs), np.max(Zabs) plt.figure(figsize = (15, 7)) plt.jet() plt.semilogy(mu1, Zabs) for l_ in approx.trainedModel.getPoles(marginal): plt.plot([np.real(l_ ** exp)] * 2, [Zmin, Zmax], 'b--') plt.plot(mu2x, [Zmax] * len(mu2x), 'kx') plt.plot([murangeExp[0][0]] * 2, [Zmin, Zmax], 'm:') plt.plot([murangeExp[1][0]] * 2, [Zmin, Zmax], 'm:') plt.xlim(mu1[0], mu1[-1]) plt.title("|Q(mu)|") plt.grid() plt.show() return mus, Z if nNone == 2: exps = [exps[m] for m in mNone] 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]]] 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)] - 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: + poles1, poles2 = [], [] + for m in mus: + m[mNone[0]] = None + pls = approx.getPoles(m) ** exps[0] + pls[np.imag(pls) > polesImTol] = np.nan + pls = np.real(pls) + poles1 += list(pls) + poles2 += [m[mNone[1]] ** exps[1]] * len(pls) + try: + 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) plt.figure(figsize = (15, 7)) plt.jet() - p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap, - levels = np.linspace(Zabsmin, Zabsmax, 50)) - if clip > 0: - plt.contour(Mu1, Mu2, Zabs, [Zabsmin]) - plt.contour(Mu1, Mu2, np.real(Z), [0.], linestyles = 'dashed') - plt.contour(Mu1, Mu2, np.imag(Z), [0.], linewidths = 1, - linestyles = 'dotted') + 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(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) + 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 if nNone == 3: exps = [exps[m] for m in mNone] - if hasattr(approx, "mus"): - mu2x = approx.mus(mNone[0]) ** exps[0] - mu2y = approx.mus(mNone[1]) ** exps[1] - mu2z = approx.mus(mNone[2]) ** exps[2] - else: - mu2x = mu0[mNone[0]] ** exps[0] - mu2y = mu0[mNone[1]] ** exps[1] - mu2z = mu0[mNone[2]] ** exps[2] murangeExp = [[murange[0][mNone[0]] ** exps[0], murange[0][mNone[1]] ** exps[1], murange[0][mNone[2]] ** exps[2]], [murange[1][mNone[0]] ** exps[0], murange[1][mNone[1]] ** exps[1], murange[1][mNone[2]] ** exps[2]]] 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) mu3 = np.linspace(murangeEff[0][mNone[2]] ** exps[2], murangeEff[1][mNone[2]] ** exps[2], nSamples) mu1s = np.power(mu1, 1. / exps[0]) mu2s = np.power(mu2, 1. / exps[1]) mu3s = np.power(mu3, 1. / exps[2]) Mu1, Mu2, Mu3 = np.meshgrid(np.real(mu1), np.real(mu2), np.real(mu3), indexing = 'ij') mus = [] mBase = copy(marginal) - for m1 in mu1s: - mBase[mNone[0]] = m1 - for m2 in mu2s: - mBase[mNone[1]] = m2 - for m3 in mu3s: - mBase[mNone[2]] = m3 - mus += [copy(mBase)] - Z = approx.trainedModel.getQVal(mus).reshape(Mu1.shape) + for m2 in mu2s: + mBase[mNone[1]] = m2 + for m3 in mu3s: + mBase[mNone[2]] = m3 + mus += [copy(mBase)] - pts = [] - for jdir in range(len(Mu3)): - z = Mu3[0, 0, jdir] - xis, yis, fis = Mu1[:, :, jdir], Mu2[:, :, jdir], Z[:, :, jdir] - plt.figure() - CS = plt.contour(xis, yis, np.real(fis), levels = [0.]) - plt.close() - for i, CSc in enumerate(CS.allsegs): - if np.isclose(CS.cvalues[i], 0.): - for j, CScj in enumerate(CSc): - for k in range(len(CScj)): - pts += [[CScj[k, 0], CScj[k, 1], z]] - pts += [[np.nan] * 3] - pts = np.array(pts) + poles1, poles2, poles3 = [], [], [] + for m in mus: + m[mNone[0]] = None + pls = approx.getPoles(m) ** exps[0] + pls[np.imag(pls) > polesImTol] = np.nan + pls = np.real(pls) + pls[pls < murangeExp[0][0]] = np.nan + pls[pls > murangeExp[1][0]] = np.nan + poles1 += list(pls) + poles2 += [m[mNone[1]] ** exps[1]] * len(pls) + poles3 += [m[mNone[2]] ** exps[2]] * len(pls) + pts = np.empty((len(poles1), 3)) + pts[:, 0] = poles1 + pts[:, 1] = poles2 + pts[:, 2] = poles3 if np.size(pts) > 0: from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize = (8, 6)) ax = Axes3D(fig) - ax.plot(pts[:, 0], pts[:, 1], pts[:, 2], 'k-') - ax.scatter(np.real(mu2x), np.real(mu2y), np.real(mu2z), '.') + ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], '.', s = 1) ax.set_xlim(murangeExp[0][0], murangeExp[1][0]) ax.set_ylim(murangeExp[0][1], murangeExp[1][1]) ax.set_zlim(murangeExp[0][2], murangeExp[1][2]) plt.show() plt.close() return pts diff --git a/examples/greedy/squareBubbleHomog.py b/examples/greedy/squareBubbleHomog.py index dbe5516..c639459 100644 --- a/examples/greedy/squareBubbleHomog.py +++ b/examples/greedy/squareBubbleHomog.py @@ -1,117 +1,117 @@ import numpy as np from rrompy.hfengines.linear_problem import \ HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.distributed_greedy import \ RationalInterpolantGreedy as Pade from rrompy.reduction_methods.distributed_greedy import \ RBDistributedGreedy as RB -from rrompy.utilities.base import squareResonances +from rrompy.utilities.numerical import squareResonances from rrompy.solver.fenics import L2NormMatrix verb = 2 timed = False algo = "Pade" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" errorEstimatorKind = "BARE" #errorEstimatorKind = "BASIC" #errorEstimatorKind = "EXACT" k0s = np.power(np.linspace(95, 149, 250), .5) #k0s = np.power(np.linspace(95, 129, 100), .5) #k0s = np.power(np.linspace(95, 109, 100), .5) k0 = np.mean(np.power(k0s, 2.)) ** .5 kl, kr = min(k0s), max(k0s) polesexact = np.unique(np.power(squareResonances(kl**2., kr**2., False), .5)) params = {'muBounds':[kl, kr], 'nTestPoints':500, 'Delta':0, 'greedyTol':1e-2, 'S':10, 'polybasis':polyBasis, 'errorEstimatorKind':errorEstimatorKind, 'interactive':False} if timed: verb = 0 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, verbosity = verb) solver.omega = np.real(k0) if algo == "Pade": approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) else: params.pop('Delta') params.pop('polybasis') params.pop('errorEstimatorKind') approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.initEstimatorNormEngine(L2NormMatrix(solver.V)) if timed: from time import clock start_time = clock() approx.greedy() print("Time: ", clock() - start_time) else: approx.greedy(True) approx.samplingEngine.verbosity = 0 approx.trainedModel.verbosity = 0 approx.verbosity = 0 from matplotlib import pyplot as plt normApp = np.zeros_like(k0s) norm = np.zeros_like(k0s) res = np.zeros_like(k0s) err = np.zeros_like(k0s) for j in range(len(k0s)): normApp[j] = approx.normApprox(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j])) / approx.estimatorNormEngine.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.semilogy(k0s, norm) plt.semilogy(k0s, normApp, '--') plt.semilogy(polesexact, 2.*np.max(norm)*np.ones_like(polesexact, dtype = float), 'm.') plt.semilogy(np.real(approx.mus(0)), 4.*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(polesexact, 2.*np.max(resApp)*np.ones_like(polesexact, dtype = float), 'm.') 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.semilogy(polesexact, 2.*np.max(err)*np.ones_like(polesexact, dtype = float), 'm.') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.plot(np.real(polesexact), np.imag(polesexact), 'm.') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/pod/PolesCentered.py b/examples/pod/PolesCentered.py index 393e6a9..ab2ee75 100644 --- a/examples/pod/PolesCentered.py +++ b/examples/pod/PolesCentered.py @@ -1,74 +1,74 @@ import numpy as np from rrompy.hfengines.linear_problem import \ HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.distributed import RationalInterpolant as RI from rrompy.reduction_methods.distributed import RBDistributed as RB -from rrompy.utilities.base import squareResonances +from rrompy.utilities.numerical import squareResonances from rrompy.parameter.parameter_sampling import ManualSampler as MS verb = 0 k0 = (12+0.j) ** .5 krange = [[9. ** .5], [15. ** .5]] Nmin, Nmax = 2, 10 Nvals = np.arange(Nmin, Nmax + 1, 2) rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'N':Nmin, 'M':0, 'S':Nmin + 1, 'POD':True, 'sampler': MS(krange, points = [k0], scaling = rescaling, scalingInv = rescalingInv)} #boolCon = lambda x : np.abs(np.imag(x)) < 1e-1 * np.abs(np.real(x) # - np.real(z0)) #cleanupParameters = {'boolCondition':boolCon, 'residueCheck':True} solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 25, verbosity = verb) solver.omega = np.real(k0) approxP = RI(solver, mu0 = k0, approxParameters = params, verbosity = verb) approxR = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) rP, rE = [None] * len(Nvals), [None] * len(Nvals) verbose = 1 for j, N in enumerate(Nvals): if verbose > 0: print('N = E = {}'.format(N)) approxP.approxParameters = {'N':N, 'S':N+1} approxR.approxParameters = {'R':N, 'S':N+1} if verbose > 1: print(approxP.approxParameters) print(approxR.approxParameters) rP[j] = approxP.getPoles() rE[j] = approxR.getPoles() if verbose > 2: print(rP) print(rE) from matplotlib import pyplot as plt plotRows = int(np.ceil(len(Nvals) / 3)) fig, axes = plt.subplots(plotRows, 3, figsize = (15, 3.5 * plotRows)) for j, N in enumerate(Nvals): i1, i2 = int(np.floor(j / 3)), j % 3 axes[i1, i2].set_title('N = E = {}'.format(N)) axes[i1, i2].plot(np.real(rP[j]), np.imag(rP[j]), 'Xb', label="Pade'", markersize = 8) axes[i1, i2].plot(np.real(rE[j]), np.imag(rE[j]), 'Pr', label="RB", markersize = 8) axes[i1, i2].axhline(linewidth=1, color='k') xmin, xmax = axes[i1, i2].get_xlim() height = (xmax - xmin) / 2. res = np.power(squareResonances(xmin**2., xmax**2., False), .5) axes[i1, i2].plot(res, np.zeros_like(res), 'ok', markersize = 4) axes[i1, i2].plot(np.real(k0), np.imag(k0), 'om', markersize = 5) axes[i1, i2].plot(np.real(k0) * np.ones(2), 1.5 * height * np.arange(-1, 3, 2), '--m') axes[i1, i2].grid() axes[i1, i2].set_xlim(xmin, xmax) axes[i1, i2].set_ylim(- height, height) p = axes[i1, i2].legend() plt.tight_layout() for j in range((len(Nvals) - 1) % 3 + 1, 3): axes[plotRows - 1, j].axis('off') diff --git a/examples/pod/PolesDistributed.py b/examples/pod/PolesDistributed.py index ee0987b..3298b77 100644 --- a/examples/pod/PolesDistributed.py +++ b/examples/pod/PolesDistributed.py @@ -1,48 +1,48 @@ from matplotlib import pyplot as plt import numpy as np from rrompy.hfengines.linear_problem import \ HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.reduction_methods.distributed import RationalInterpolant as RI from rrompy.parameter.parameter_sampling import QuadratureSampler as QS -from rrompy.utilities.base import squareResonances +from rrompy.utilities.numerical import squareResonances verb = 0 ks = [1., 46. ** .5] solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, verbosity = verb) k0 = np.mean(np.power(ks, 2.)) ** .5 k0 = 3.46104724 solver.omega = np.real(k0) rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) nsets = 12 paramsPade = {'S':2, 'POD':True, 'polybasis':"LEGENDRE", 'sampler':QS(ks, "UNIFORM", rescaling, rescalingInv)} approx = RI(solver, mu0 = k0, approxParameters = paramsPade, verbosity = verb) poles = [None] * nsets samples = [None] * nsets polesexact = np.unique(np.power( squareResonances(ks[0]**2., ks[1]**2., False), .5)) for i in range(1, nsets + 1): print("N = {}".format(4 * i)) approx.approxParameters = {'N': 4 * i, 'M': 4 * i, 'S': 4 * i + 1} approx.setupApprox() poles[i - 1] = approx.getPoles() samples[i - 1] = (approx.mus ** 2.).data for i in range(1, nsets + 1): plt.figure() plt.plot(np.real(poles[i - 1]), np.imag(poles[i - 1]), 'kx') plt.plot(polesexact, np.zeros_like(polesexact), 'm.') plt.plot(np.real(samples[i - 1]), np.imag(samples[i - 1]), 'r*') plt.xlim(ks) plt.ylim((ks[0] - ks[1]) / 2., (ks[1] - ks[0]) / 2.) plt.title("N = {}, Neff = {}".format(4 * i, len(poles[i - 1]))) plt.grid() plt.show() plt.close()