diff --git a/examples/2d/greedy/fracture_pod.py b/examples/2d/greedy/fracture_pod.py index f752364..55b2181 100644 --- a/examples/2d/greedy/fracture_pod.py +++ b/examples/2d/greedy/fracture_pod.py @@ -1,220 +1,220 @@ import numpy as np from membrane_fracture_engine 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 = 16 show_sample = False show_norm = True MN = 1 S = (MN + 2) * (MN + 1) // 2 greedyTol = 1e-1 collinearityTol = 0. nTestPoints = 900 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 = 'EIM_DIAGONAL' errorEstimatorKind = 'EIM_INTERPOLATORY' 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.greedy(True) +approx.setupApprox(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 from plot_zero_set import plotZeroSet2 muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0, 200, [2., 1.]) if show_norm: from plot_inf_set import plotInfSet2 muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2( murange, murangeEff, approx, mu0, 50, [2., 1.], relative = True, nobeta = True) try: QV = approx.trainedModel.getQVal(muInfVals) 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) * QV).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() except: pass approx.trainedModel.verbosity = verb try: print(np.sort(approx.getPoles([None, .5]) ** 2.)) except: pass diff --git a/examples/airfoil/greedy.py b/examples/airfoil/greedy.py index 2f38d0a..16cd66c 100644 --- a/examples/airfoil/greedy.py +++ b/examples/airfoil/greedy.py @@ -1,103 +1,103 @@ import numpy as np from airfoil_engine import AirfoilScatteringEngine 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 verb = 5 timed = False algo = "RI" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" k0s = np.linspace(5, 10, 100) k0 = np.mean(k0s) kl, kr = min(k0s), max(k0s) params = {'S':5, 'sampler':QS([kl, kr], "UNIFORM"), 'nTestPoints':500, 'greedyTol':1e-2, 'polybasis':polyBasis, 'errorEstimatorKind':'INTERPOLATORY'} ######### kappa = 10 theta = np.pi * - 45 / 180. solver = AirfoilScatteringEngine(kappa, theta, verbosity = verb, degree_threshold = 8) ######### 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() + approx.setupApprox() print("Time: ", clock() - start_time) else: - approx.greedy(True) + approx.setupApprox(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 kl, kr = np.real(kl), np.real(kr) from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApprox(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j])) / approx.estimatorNormEngine.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.plot(k0s, norm) plt.plot(k0s, normApp, '--') plt.plot(np.real(approx.mus.data), 1.05*np.max(norm)*np.ones_like(approx.mus.data, 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.data), 4.*np.max(resApp)*np.ones_like(approx.mus.data, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/all_forcing/greedy.py b/examples/all_forcing/greedy.py index 94e5446..e54d8c9 100644 --- a/examples/all_forcing/greedy.py +++ b/examples/all_forcing/greedy.py @@ -1,94 +1,94 @@ import numpy as np from all_forcing_engine import AllForcingEngine 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 verb = 5 timed = False algo = "RI" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" if timed: verb = 0 z0s = np.linspace(-3., 3., 100) z0 = np.mean(z0s) zl, zr = min(z0s), max(z0s) n = 30 solver = AllForcingEngine(mu0 = z0, n = n, degree_threshold = 8, verbosity = 0) params = {'sampler':QS([zl, zr], "UNIFORM"), 'nTestPoints':500, 'greedyTol':1e-2, 'S':2, 'polybasis':polyBasis, # 'errorEstimatorKind':'DIAGONAL'} # 'errorEstimatorKind':'INTERPOLATORY'} 'errorEstimatorKind':'AFFINE'} if algo == "RI": approx = RI(solver, mu0 = solver.mu0, approxParameters = params, verbosity = verb) else: params.pop("polybasis") params.pop("errorEstimatorKind") approx = RB(solver, mu0 = solver.mu0, approxParameters = params, verbosity = verb) approx.initEstimatorNormEngine(L2NormMatrix(solver.V)) if timed: from time import clock start_time = clock() - approx.greedy() + approx.setupApprox() print("Time: ", clock() - start_time) else: - approx.greedy(True) + approx.setupApprox(True) polesApp = approx.getPoles() print("Poles:\n", polesApp) approx.samplingEngine.verbosity = 0 approx.trainedModel.verbosity = 0 approx.verbosity = 0 zl, zr = np.real(zl), np.real(zr) from matplotlib import pyplot as plt normApp = np.zeros(len(z0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(z0s)): normApp[j] = approx.normApprox(z0s[j]) norm[j] = approx.normHF(z0s[j]) res[j] = (approx.estimatorNormEngine.norm(approx.getRes(z0s[j])) / approx.estimatorNormEngine.norm(approx.getRHS(z0s[j]))) err[j] = approx.normErr(z0s[j]) / norm[j] resApp = approx.errorEstimator(z0s) plt.figure() plt.semilogy(z0s, norm) plt.semilogy(z0s, normApp, '--') plt.semilogy(np.real(approx.mus.data), 1.05*np.max(norm)*np.ones_like(approx.mus.data, dtype = float), 'rx') plt.xlim([zl, zr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(z0s, res) plt.semilogy(z0s, resApp, '--') plt.semilogy(np.real(approx.mus.data), approx.greedyTol*np.ones_like(approx.mus.data, dtype = float), 'rx') plt.xlim([zl, zr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(z0s, err) plt.xlim([zl, zr]) plt.grid() plt.show() plt.close() diff --git a/examples/diapason/greedy.py b/examples/diapason/greedy.py index cc85744..35cfcfb 100644 --- a/examples/diapason/greedy.py +++ b/examples/diapason/greedy.py @@ -1,127 +1,127 @@ import numpy as np from diapason_engine import DiapasonEngine, DiapasonEngineDamped 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 verb = 5 timed = False algo = "RI" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" if timed: verb = 0 dampingEta = 0 * 1e4 / 2. / np.pi k0s = np.linspace(2.5e2, 7.5e3, 100) #k0s = np.linspace(2.5e3, 1.5e4, 100) #k0s = np.linspace(5.0e4, 1.0e5, 100) k0s = np.linspace(2.0e5, 2.5e5, 100) kl, kr = min(k0s), max(k0s) theta = 20. * np.pi / 180. phi = 10. * np.pi / 180. c = 3.e2 rho = 8e3 * (2. * np.pi) ** 2. E = 1.93e11 nu = .3 T = 1e6 ### if np.isclose(dampingEta, 0.): solver = DiapasonEngine(kappa = np.mean(np.power(k0s, 2.)) ** .5, c = c, rho = rho, E = E, nu = nu, T = T, theta = theta, phi = phi, meshNo = 1, degree_threshold = 8, verbosity = 0) else: solver = DiapasonEngineDamped(kappa = np.mean(k0s), c = c, rho = rho, E = E, nu = nu, T = T, theta = theta, phi = phi, dampingEta = dampingEta, meshNo = 1, degree_threshold = 8, verbosity = 0) params = {'sampler':QS([kl, kr], "UNIFORM"),#, solver.rescalingExp), 'nTestPoints':500, 'greedyTol':1e-2, 'S':5, 'polybasis':polyBasis, # 'errorEstimatorKind':'DIAGONAL'} # 'errorEstimatorKind':'INTERPOLATORY'} 'errorEstimatorKind':'AFFINE'} if algo == "RI": approx = RI(solver, mu0 = solver.mu0, approxParameters = params, verbosity = verb) else: params.pop("polybasis") params.pop("errorEstimatorKind") approx = RB(solver, mu0 = solver.mu0, approxParameters = params, verbosity = verb) approx.initEstimatorNormEngine(L2NormMatrix(solver.V)) if timed: from time import clock start_time = clock() - approx.greedy() + approx.setupApprox() print("Time: ", clock() - start_time) else: - approx.greedy(True) + approx.setupApprox(True) polesApp = approx.getPoles() print("Poles:\n", polesApp) approx.samplingEngine.verbosity = 0 approx.trainedModel.verbosity = 0 approx.verbosity = 0 kl, kr = np.real(kl), np.real(kr) from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApprox(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j])) / approx.estimatorNormEngine.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / norm[j] resApp = approx.errorEstimator(k0s) plt.figure() plt.semilogy(k0s, norm) plt.semilogy(k0s, normApp, '--') plt.semilogy(np.real(approx.mus.data), 1.05*np.max(norm)*np.ones_like(approx.mus.data, 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.data), approx.greedyTol*np.ones_like(approx.mus.data, dtype = float), 'rx') plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() plt.figure() plt.semilogy(k0s, err) plt.xlim([kl, kr]) plt.grid() plt.show() plt.close() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesAppEff = polesApp[~mask] plt.figure() plt.plot(np.real(polesAppEff), np.imag(polesAppEff), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/from_papers/greedy_internalBox.py b/examples/from_papers/greedy_internalBox.py index 26e10a1..5dd6fbe 100644 --- a/examples/from_papers/greedy_internalBox.py +++ b/examples/from_papers/greedy_internalBox.py @@ -1,122 +1,122 @@ 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':1e-2, 'S':2, 'polybasis':polyBasis, 'collinearityTol':1e3, # 'errorEstimatorKind':'AFFINE'} # 'errorEstimatorKind':'DISCREPANCY'} # 'errorEstimatorKind':'INTERPOLATORY'} # 'errorEstimatorKind':'NONE'} 'errorEstimatorKind':'LOOK_AHEAD'} 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() + approx.setupApprox() print("Time: ", clock() - start_time) else: - approx.greedy(True) + approx.setupApprox(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 kl, kr = np.real(kl), np.real(kr) from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApprox(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j])) / approx.estimatorNormEngine.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.plot(k0s, norm) plt.plot(k0s, normApp, '--') plt.plot(np.real(approx.mus(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) if algo == "RB" or params['errorEstimatorKind'] != 'LOOK_AHEAD': 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) if algo == "RI" and params['errorEstimatorKind'] == 'LOOK_AHEAD': 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() polesApp = approx.getPoles() mask = (np.real(polesApp) < kl) | (np.real(polesApp) > kr) print("Outliers:", polesApp[mask]) polesApp = polesApp[~mask] plt.figure() plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') plt.axis('equal') plt.grid() plt.show() plt.close() diff --git a/examples/from_papers/greedy_scatteringAirfoil.py b/examples/from_papers/greedy_scatteringAirfoil.py index 00dafc2..0e0db20 100644 --- a/examples/from_papers/greedy_scatteringAirfoil.py +++ b/examples/from_papers/greedy_scatteringAirfoil.py @@ -1,147 +1,147 @@ import numpy as np import fenics as fen import ufl from helmholtz_box_scattering_problem_engine import \ HelmholtzBoxScatteringProblemEngine as HSP from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RI from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RB from rrompy.solver.fenics import fenONE, L2NormMatrix from rrompy.parameter.parameter_sampling import QuadratureSampler as QS verb = 2 timed = False algo = "RI" #algo = "RB" polyBasis = "LEGENDRE" #polyBasis = "CHEBYSHEV" #polyBasis = "MONOMIAL" k0s = np.linspace(5, 20, 25) k0 = np.mean(k0s) kl, kr = min(k0s), max(k0s) params = {'sampler':QS([kl, kr], "UNIFORM"), 'nTestPoints':500, 'greedyTol':1e-2, 'S':10, 'polybasis':polyBasis, 'errorEstimatorKind':'DIAGONAL'} # 'errorEstimatorKind':'INTERPOLATORY'} # 'errorEstimatorKind':'AFFINE'} ######### PI = np.pi R = 2 def Dboundary(x, on_boundary): return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R kappa = 10 theta = PI * - 45 / 180. mu = 1.1 epsilon = .1 mesh = fen.Mesh('../data/mesh/airfoil.xml') c, s = np.cos(theta), np.sin(theta) x, y = fen.SpatialCoordinate(mesh)[:] u0R = - fen.cos(kappa * (c * x + s * y)) u0I = - fen.sin(kappa * (c * x + s * y)) checkReal = x**2-x+y**2 rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25 phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2 phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2 kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/ ((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.) )**.5 - mu kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/ ((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.) )**.5 - mu Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1 Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1 cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE) c_F = fen.Constant(.1) cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE) c_F = fen.Constant(.1) cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F) cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F) a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF) solver = HSP(R, np.real(k0), theta, n = 1, verbosity = verb, degree_threshold = 8) solver.omega = np.real(k0) solver.V = fen.FunctionSpace(mesh, "P", 2) solver.diffusivity = a solver.DirichletBoundary = Dboundary solver.RobinBoundary = "REST" solver.DirichletDatum = [u0R, u0I] ######### 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() + approx.setupApprox() print("Time: ", clock() - start_time) else: - approx.greedy(True) + approx.setupApprox(True) approx.samplingEngine.verbosity = 0 approx.verbosity = 0 kl, kr = np.real(kl), np.real(kr) from matplotlib import pyplot as plt normApp = np.zeros(len(k0s)) norm = np.zeros_like(normApp) res = np.zeros_like(normApp) err = np.zeros_like(normApp) for j in range(len(k0s)): normApp[j] = approx.normApprox(k0s[j]) norm[j] = approx.normHF(k0s[j]) res[j] = (approx.estimatorNormEngine.norm(approx.getRes(k0s[j])) / approx.estimatorNormEngine.norm(approx.getRHS(k0s[j]))) err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) resApp = approx.errorEstimator(k0s) plt.figure() plt.plot(k0s, norm) plt.plot(k0s, normApp, '--') plt.plot(np.real(approx.mus(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()