Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/examples/2d/greedy/fracture_pod.py
3ed8cef19cfcmaster

Last Change
Raw File

fracture_pod.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem.bidimensional import \
MembraneFractureEngine as MFE
from rrompy.reduction_methods.greedy import RationalInterpolantGreedy as RIG
from rrompy.reduction_methods.greedy import ReducedBasisGreedy as RBG
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST)
verb = 15
size = 2
show_sample = True
show_norm = True
MN = 1
S = (MN + 2) * (MN + 1) // 2
greedyTol = 1e-2
collinearityTol = 0.
nTestPoints = 400
algo = "rational"
algo = "RB"
if algo == "rational":
radial = ""
#radial = "_gaussian"
#radial = "_thinplate"
#radial = "_multiquadric"
rW0 = 5.
radialWeight = [rW0] * 2
polybasis = "CHEBYSHEV"
polybasis = "LEGENDRE"
#polybasis = "MONOMIAL"
errorEstimatorKind = 'AFFINE'
errorEstimatorKind = 'DISCREPANCY'
# errorEstimatorKind = 'INTERPOLATORY'
# errorEstimatorKind = 'LOCALL2'
# errorEstimatorKind = 'DIAGONAL'
if size == 1: # below
mu0 = [40 ** .5, .4]
mutar = [45 ** .5, .4]
murange = [[30 ** .5, .3], [50 ** .5, .5]]
elif size == 2: # top
mu0 = [40 ** .5, .6]
mutar = [45 ** .5, .6]
murange = [[30 ** .5, .5], [50 ** .5, .7]]
elif size == 3: # interesting
mu0 = [40 ** .5, .5]
mutar = [45 ** .5, .5]
murange = [[30 ** .5, .3], [50 ** .5, .7]]
elif size == 4: # wide_low
mu0 = [40 ** .5, .2]
mutar = [45 ** .5, .2]
murange = [[10 ** .5, .1], [70 ** .5, .3]]
elif size == 5: # wide_hi
mu0 = [40 ** .5, .8]
mutar = [45 ** .5, .8]
murange = [[10 ** .5, .7], [70 ** .5, .9]]
elif size == 6: # top_zoom
mu0 = [50 ** .5, .8]
mutar = [55 ** .5, .8]
murange = [[40 ** .5, .7], [60 ** .5, .9]]
elif size == 7: # huge
mu0 = [50 ** .5, .5]
mutar = [55 ** .5, .5]
murange = [[10 ** .5, .2], [90 ** .5, .8]]
elif size == 11: #L below
mu0 = [110 ** .5, .4]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .3], [130 ** .5, .5]]
elif size == 12: #L top
mu0 = [110 ** .5, .6]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .7]]
elif size == 13: #L interesting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .5]
murange = [[90 ** .5, .3], [130 ** .5, .7]]
elif size == 14: #L belowbelow
mu0 = [110 ** .5, .2]
mutar = [115 ** .5, .2]
murange = [[90 ** .5, .1], [130 ** .5, .3]]
elif size == 15: #L toptop
mu0 = [110 ** .5, .8]
mutar = [115 ** .5, .8]
murange = [[90 ** .5, .7], [130 ** .5, .9]]
elif size == 16: #L interestinginteresting
mu0 = [110 ** .5, .5]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .1], [130 ** .5, .9]]
elif size == 17: #L interestingtop
mu0 = [110 ** .5, .7]
mutar = [115 ** .5, .6]
murange = [[90 ** .5, .5], [130 ** .5, .9]]
elif size == 18: #L interestingbelow
mu0 = [110 ** .5, .3]
mutar = [115 ** .5, .4]
murange = [[90 ** .5, .1], [130 ** .5, .5]]
elif size == 100: # tiny
mu0 = [32.5 ** .5, .5]
mutar = [34 ** .5, .5]
murange = [[30 ** .5, .3], [35 ** .5, .7]]
aEff = 1#.05
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
aEff*murange[0][1] + bEff*murange[1][1]],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
aEff*murange[1][1] + bEff*murange[0][1]]]
H = 1.
L = .75
delta = .05
n = 20
solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
rescalingExp = [2., 1.]
if algo == "rational":
params = {'S':S, 'POD':True, 'greedyTol':greedyTol,
'sampler':QS(murange, 'UNIFORM', scalingExp = rescalingExp),
'nTestPoints':nTestPoints, 'polybasis':polybasis + radial,
'radialDirectionalWeights':radialWeight,
'errorEstimatorKind':errorEstimatorKind,
'trainSetGenerator':QST(murange, 'CHEBYSHEV',
scalingExp = rescalingExp)}
method = RIG
else: #if algo == "RB":
params = {'S':S, 'POD':True, 'greedyTol':greedyTol,
'sampler':QS(murange, 'UNIFORM', scalingExp = rescalingExp),
'nTestPoints':nTestPoints,
'trainSetGenerator':QST(murange, 'CHEBYSHEV',
scalingExp = rescalingExp)}
method = RBG
approx = method(solver, mu0 = mu0, approxParameters = params, verbosity = verb)
approx.samplingEngine.allowRepeatedSamples = False
approx.greedy(True)
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = mutar[1]
y = fen.SpatialCoordinate(solver.V.mesh())[1]
warp1, warpI1 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. * L]]))
warp2, warpI2 = affine_warping(solver.V.mesh(),
np.array([[1, 0], [0, 2. - 2. * L]]))
warp = ufl.conditional(ufl.ge(y, 0.), warp1, warp2)
warpI = ufl.conditional(ufl.ge(y, 0.), warpI1, warpI2)
approx.plotApprox(mutar, [warp, warpI], name = 'u_app', what = "REAL")
approx.plotHF(mutar, [warp, warpI], name = 'u_HF', what = "REAL")
approx.plotErr(mutar, [warp, warpI], name = 'err', what = "REAL")
# approx.plotRes(mutar, [warp, warpI], name = 'res', what = "REAL")
appErr = approx.normErr(mutar)
solNorm = approx.normHF(mutar)
resNorm = approx.normRes(mutar)
RHSNorm = approx.normRHS(mutar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
verb = approx.trainedModel.verbosity
approx.trainedModel.verbosity = 5
try:
from plot_zero_set import plotZeroSet2
muZeroVals, Qvals = plotZeroSet2(murange, murangeEff, approx, mu0,
200, [2., 1.])
except: pass
if show_norm:
solver._solveBatchSize = 100
from plot_inf_set import plotInfSet2
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet2(
murange, murangeEff, approx, mu0, 50,
[2., 1.], relative = True,
nobeta = True)
import warnings
from matplotlib import pyplot as plt
mu2x, mu2y = approx.mus(0) ** 2., approx.mus(1) ** 1.
murangeExp = [[murange[0][0] ** 2., murange[0][1]],
[murange[1][0] ** 2., murange[1][1]]]
mu1s = np.unique([m[0] for m in muInfVals])
mu2s = np.unique([m[1] for m in muInfVals])
mu1 = np.power(mu1s, 2.)
mu2 = np.power(mu2s, 1.)
Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
Reste = (approx.errorEstimator(muInfVals)
* approx.trainedModel.getQVal(muInfVals)).reshape(normEx.shape)
Rest = np.log10(Reste)
Restmin, Restmax = np.min(Rest), np.max(Rest)
cmap = plt.cm.jet
warnings.simplefilter("ignore", category = (UserWarning,
np.ComplexWarning))
plt.figure(figsize = (15, 7))
plt.jet()
p = plt.contourf(Mu1, Mu2, Rest,
levels = np.linspace(Restmin, Restmax, 50))
plt.plot(np.real(mu2x), np.real(mu2y), 'kx')
for j, xy in enumerate(zip(np.real(mu2x), np.real(mu2y))):
plt.annotate("{}".format(j), xy)
plt.plot([murangeExp[0][0]] * 2,
[murangeExp[0][1], murangeExp[1][1]], 'm:')
plt.plot([murangeExp[0][0], murangeExp[1][0]],
[murangeExp[1][1]] * 2, 'm:')
plt.plot([murangeExp[1][0]] * 2,
[murangeExp[1][1], murangeExp[0][1]], 'm:')
plt.plot([murangeExp[1][0], murangeExp[0][0]],
[murangeExp[0][1]] * 2, 'm:')
plt.title("log10|res_est(mu)|")
plt.colorbar(p)
plt.grid()
plt.show()
approx.trainedModel.verbosity = verb
try:
print(np.sort(approx.getPoles([None, .5]) ** 2.))
except: pass
c4science · Help