diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py index cbc4c94..0cb10f1 100644 --- a/examples/3d/pod/fracture3_pod.py +++ b/examples/3d/pod/fracture3_pod.py @@ -1,196 +1,204 @@ 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.centered import RationalPade as RP from rrompy.reduction_methods.distributed import RationalInterpolant as RI from rrompy.reduction_methods.centered import RBCentered as RBC from rrompy.reduction_methods.distributed import RBDistributed as RBD from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, QuadratureSamplerTotal as QST, ManualSampler as MS, RandomSampler as RS) verb = 50 -size = 2 +size = 4 show_sample = True show_norm = False clip = -1 #clip = .4 #clip = .6 homogeneize = False #homogeneize = True Delta = 0 -MN = 15 +MN = 10 R = (MN + 3) * (MN + 2) * (MN + 1) // 6 S = [int(np.ceil(R ** (1. / 3.)))] * 3 PODTol = 1e-8 samples = "centered" samples = "centered_fake" samples = "distributed" algo = "rational" #algo = "RB" sampling = "quadrature" sampling = "quadrature_total" -#sampling = "random" +sampling = "random" if samples == "distributed": radial = 0 # radial = "gaussian" # radial = "thinplate" # radial = "multiquadric" rW0 = 5. radialWeight = [rW0] * 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 = 20 +n = 40 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_fake": params['polybasis'] = "MONOMIAL" params['S'] = R method = RI else: params['S'] = R method = RP 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_fake": params['S'] = R method = RBD else: params['S'] = R method = RBC if samples == "distributed": if sampling == "quadrature": params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) # params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling, # scalingInv = rescalingInv) # params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling, # scalingInv = rescalingInv) params['S'] = [max(j, MN + 1) for j in params['S']] elif sampling == "quadrature_total": params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R else: # if sampling == "random": params['sampler'] = RS(murange, "HALTON", scaling = rescaling, scalingInv = rescalingInv) params['S'] = R elif samples == "centered_fake": params['sampler'] = MS(murange, points = [mu0], scaling = rescaling, scalingInv = rescalingInv) approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb, homogeneized = homogeneize) if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False approx.setupApprox() if show_sample: 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) 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), '.') 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) - plotZeroSet3(murange, murangeEff, approx, mu0, 100, [None, None, None], - [2., 1., 1.], clip = clip) + plotZeroSet3(murange, murangeEff, approx, mu0, 50, [None, None, None], + [2., 1., 1.], clip = clip, fourD = False) 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 02e3cd2..1b3ede3 100644 --- a/examples/3d/pod/plot_zero_set_3.py +++ b/examples/3d/pod/plot_zero_set_3.py @@ -1,155 +1,167 @@ 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, imagTol = 1e-8, fourD = True): 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: cmap = plt.cm.jet warnings.simplefilter("ignore", category = UserWarning) plt.figure(figsize = (15, 7)) plt.jet() p = plt.contourf(Mu1, Mu2, Zabs, cmap = cmap, levels = np.linspace(Zabsmin, Zabsmax, 50)) if clip > 0: plt.contour(Mu1, Mu2, Zabs, [Zabsmin]) plt.contour(Mu1, Mu2, np.real(Z), [0.], linestyles = 'dashed') plt.contour(Mu1, Mu2, np.imag(Z), [0.], linewidths = 1, linestyles = 'dotted') plt.plot(mu2x, mu2y, 'kx') plt.plot([murangeExp[0][0]] * 2, [murangeExp[0][1], murangeExp[1][1]], 'm:') plt.plot([murangeExp[0][0], murangeExp[1][0]], [murangeExp[1][1]] * 2, 'm:') plt.plot([murangeExp[1][0]] * 2, [murangeExp[1][1], murangeExp[0][1]], 'm:') plt.plot([murangeExp[1][0], murangeExp[0][0]], [murangeExp[0][1]] * 2, 'm:') plt.colorbar(p) plt.title("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]) - vals = np.empty((nSamples * nSamples * approx.N, 3), - dtype = np.complex) - for i2, m2 in enumerate(mu2s): - for i3, m3 in enumerate(mu3s): - zs = approx.getPoles([None, m2, m3]) - idx = nSamples * (nSamples * np.arange(approx.N) + i2) + i3 - vals[idx, :] = [[z, m2, m3] for z in zs ** exps[0]] - scaled_imag = np.imag(vals[:, 0]) / imagTol - rgba_colors = np.zeros((len(vals), 4)) - rgba_colors[:, 0] = np.arctan(.1 * scaled_imag) / np.pi + .5 - rgba_colors[:, 1] = - np.arctan(.1 * scaled_imag) / np.pi + .5 - rgba_colors[:, 2] = rgba_colors[:, 0] - rgba_colors[:, 3] = (np.abs(scaled_imag) ** 2. + 1) ** -.5 - vals = np.real(vals) - vals[vals[:, 0] < murangeExp[0][0], 0] = np.nan - vals[vals[:, 0] > murangeExp[1][0], 0] = np.nan + 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) + + 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) from mpl_toolkits.mplot3d import Axes3D - warnings.simplefilter("ignore", category = UserWarning) fig = plt.figure(figsize = (8, 6)) - plt.jet() ax = Axes3D(fig) - ax.scatter(vals[:, 0], vals[:, 1], vals[:, 2], - marker = 'o', s = 5, c = rgba_colors) - ax.scatter(mu2x, mu2y, mu2z, 'k.') + ax.plot(pts[:, 0], pts[:, 1], pts[:, 2], 'k-') + ax.scatter(mu2x, mu2y, mu2z, '.') + plt.show() + plt.close() 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() -