Page MenuHomec4science
Files F60503243

fracture_pod_nodomain.py
No OneTemporary
Actions

File Metadata

Created
Tue, Apr 30, 17:02

fracture_pod_nodomain.py
View Options

import numpy as np
from rrompy.hfengines.linear_problem import MembraneFractureEngineNoDomain \
as MFEND
from rrompy.reduction_methods.standard import RationalInterpolant as RI
from rrompy.reduction_methods.standard import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
QuadratureSamplerTotal as QST,
ManualSampler as MS,
RandomSampler as RS)
verb = 5
size = 1
show_sample = True
show_norm = True
ignore_forcing = True
ignore_forcing = False
clip = -1
#clip = .4
#clip = .6
homogeneize = False
#homogeneize = True
Delta = 0
MN = 6
R = MN + 1
samples = "centered"
samples = "standard"
algo = "rational"
#algo = "RB"
sampling = "quadrature"
sampling = "quadrature_total"
sampling = "random"
if samples == "standard":
radial = ""
# radial = "_gaussian"
# radial = "_thinplate"
# radial = "_multiquadric"
radialWeight = [2.]
assert Delta <= 0
if size == 1: # below
mu0Aug = [40 ** .5, .4]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 2: # top
mu0Aug = [40 ** .5, .6]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 3: # interesting
mu0Aug = [40 ** .5, .5]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[30 ** .5], [50 ** .5]]
elif size == 4: # wide_low
mu0Aug = [40 ** .5, .2]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 5: # wide_hi
mu0Aug = [40 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 45 ** .5
murange = [[10 ** .5], [70 ** .5]]
elif size == 6: # top_zoom
mu0Aug = [50 ** .5, .8]
mu0 = mu0Aug[0]
mutar = 55 ** .5
murange = [[40 ** .5], [60 ** .5]]
aEff = 1.#25
bEff = 1. - aEff
murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5],
[(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5]]
H = 1.
L = .75
delta = .05
n = 20
solver = MFEND(mu0 = mu0Aug, H = H, L = L, delta = delta, n = n,
verbosity = verb, homogeneized = homogeneize)
rescalingExp = 2.
if algo == "rational":
params = {'N':MN, 'M':MN + Delta, 'S':R, 'POD':True}
if samples == "standard":
params['polybasis'] = "CHEBYSHEV" + radial
# params['polybasis'] = "LEGENDRE" + radial
# params['polybasis'] = "MONOMIAL" + radial
params['radialDirectionalWeights'] = radialWeight
elif samples == "centered":
params['polybasis'] = "MONOMIAL"
method = RI
else: #if algo == "RB":
params = {'R':R, 'S':R, 'POD':True}
method = RB
if samples == "standard":
if sampling == "quadrature":
params['sampler'] = QS(murange, "CHEBYSHEV", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "GAUSSLEGENDRE", scalingExp = rescalingExp)
# params['sampler'] = QS(murange, "UNIFORM", scalingExp = rescalingExp)
elif sampling == "quadrature_total":
params['sampler'] = QST(murange, "CHEBYSHEV", scalingExp = rescalingExp)
else: # if sampling == "random":
params['sampler'] = RS(murange, "HALTON", scalingExp = rescalingExp)
elif samples == "centered":
params['sampler'] = MS(murange, points = [mu0], scalingExp = rescalingExp)
approx = method(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
if samples == "standard": approx.samplingEngine.allowRepeatedSamples = False
approx.setupApprox()
if show_sample:
import ufl
import fenics as fen
from rrompy.solver.fenics import affine_warping
L = solver.lFrac
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)))
if algo == "rational":
from plot_zero_set import plotZeroSet1
muZeroVals, Qvals = plotZeroSet1(murange, murangeEff, approx, mu0,
1000, 2.)
if show_norm:
solver._solveBatchSize = 10
from plot_inf_set import plotInfSet1
muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet1(
murange, murangeEff, approx, mu0,
200, 2., relative = False,
normalizeDen = True)
print(approx.getPoles() ** 2.)

Event Timeline

c4science · Help