Page MenuHomec4science
Files F63083013

damped_mass_chain.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 17, 16:11

damped_mass_chain.py
View Options

### example from Lohmann, Eid. Efficient Order Reduction of Parametric and
### Nonlinear Models by Superposition of Locally Reduced Models.
from copy import deepcopy as copy
import numpy as np
import matplotlib.pyplot as plt
from rrompy.hfengines.scipy_engines import (OscillatorProblemEngine,
AugmentedOscillatorProblemEngine)
from rrompy.reduction_methods import (NearestNeighbor as NN,
RationalInterpolant as RI,
RationalInterpolantGreedy as RIG,
RationalInterpolantPivoted as RIP,
RationalInterpolantGreedyPivoted as RIGP)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
EmptySampler as ES)
def bode(freq, mus, models):
x = np.hstack((freq.reshape(-1, 1), np.tile(np.array(mus).reshape(1, -1),
(len(freq), 1))))
fig = plt.figure(figsize = (7.5, 4.5))
ax = fig.add_subplot(1, 1, 1)
for model in models: ax.semilogx(x[:, 0],
20 * np.log10(np.abs(model(x).data.T)))
ax.set_xlim(freq[0], freq[-1])
ax.set_xlabel("frequency")
ax.set_ylabel("output")
ax.set_title(str(mus))
ax.grid()
plt.show()
return fig
N, order = 4, 1 #+ 1
if order == 1:
engine = AugmentedOscillatorProblemEngine
else:
engine = OscillatorProblemEngine
SMarginal = 1 + 1 #+ 1
M = [np.array([1., 5., 25., 125.])]
D = [np.zeros((N, N))]
D[0][0, 1], D[0][1, 2], D[0][2, 3], D[0][3, 3] = .1, .4, 1.6, 0.
D[0] = D[0] + D[0].T
K = [np.zeros((N, N))]
K[0][0, 1], K[0][1, 2], K[0][2, 3], K[0][0, 3], K[0][3, 3] = 9., 3., 1., 1., 2.
K[0] = K[0] + K[0].T
B, C = np.zeros((N, 1)), np.zeros((1, N))
B[0, 0], C[0, 3] = 27., 1.
if SMarginal > 1:
M += [np.zeros(N)]
D += [np.zeros((N, N))]
D[1][3, 3] = 1.
D[1] = D[1] + D[1].T
K += [np.zeros((N, N))]
K[1][0, 3], K[1][3, 3] = 2., 2.
K[1] = K[1] + K[1].T
solver = engine(M, D, K, B, C, order)
ss = [0., 10.]
s0, krange = np.mean(ss), [[ss[0]], [ss[-1]]]
k0, srange = [s0], copy(krange)
freq, mu = np.logspace(-2, 1, 100), []
if SMarginal > 1:
ms = [0., 1.]
m0, mrange = np.mean(ms), [[ms[0]], [ms[-1]]]
krange[0] += mrange[0]
krange[1] += mrange[1]
k0 += [m0]
mu = [.5 * (ms[1] - ms[0]) / (SMarginal - 1)]
for method in ["RI", "RI_GREEDY"]:
print("Testing {} method".format(method))
if method == "RI":
params = {'S':8, 'POD':True, 'polybasis':"CHEBYSHEV"}
if SMarginal > 1:
algo = RIP
else:
params['sampler'] = QS(srange, "CHEBYSHEV")
algo = RI
if method == "RI_GREEDY":
params = {'S':5, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2,
'errorEstimatorKind':"DISCREPANCY",
'trainSetGenerator':QS(srange, "CHEBYSHEV")}
if SMarginal > 1:
algo = RIGP
else:
params['sampler'] = QS(srange, "UNIFORM")
algo = RIG
if SMarginal > 1:
params["paramsMarginal"] = {"MMarginal": SMarginal - 1}
params['SMarginal'] = SMarginal
params['polybasisMarginal'] = "MONOMIAL"
params['radialDirectionalWeightsMarginal'] = [2. / (ms[1] - ms[0])]
params['matchingWeight'] = 1.
#params['cutOffTolerance'] = 2.
params['samplerPivot'] = QS(srange, "UNIFORM")
params['samplerMarginal'] = QS(mrange, "UNIFORM")
approx = algo([0], solver, mu0 = k0, approx_state = True,
approxParameters = params, verbosity = 5,
storeAllSamples = True)
else:
approx = algo(solver, mu0 = k0, approx_state = True,
approxParameters = params, verbosity = 5)
if "GREEDY" in method:
approx.setupApprox("LAST")
else:
approx.setupApprox()
approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 5,
approxParameters = {'S':len(approx.mus), 'POD':True,
'sampler':ES()})
if SMarginal > 1:
approxNN.setSamples(approx.storedSamplesFilenames)
approx.purgeStoredSamples()
for m in approx.musMarginal:
bode(freq, m[0], [approx.getHF, approx.getApprox,
approxNN.getApprox])
else:
approxNN.setSamples(approx.samplingEngine)
bode(freq, mu, [approx.getHF, approx.getApprox, approxNN.getApprox])
if SMarginal > 1:
bode(freq, [1.5 * ms[1]],
[approx.getHF, approx.getApprox, approxNN.getApprox])
bode(freq, [2. * ms[1]],
[approx.getHF, approx.getApprox, approxNN.getApprox])
verb = approx.verbosity
approx.verbosity = 0
mspace = np.linspace(ms[0], ms[-1], 10)
for j, t in enumerate(mspace):
pls = approx.getPoles([None, t])
if j == 0:
poles = np.empty((len(mspace), len(pls)), dtype = np.complex)
poles[j] = pls
for j, t in enumerate(approx.musMarginal):
pls = approx.getPoles([None, t[0][0]])
if j == 0:
polesE = np.empty((SMarginal, len(pls)), dtype = np.complex)
polesE[j] = pls
approx.verbosity = verb
fig = plt.figure(figsize = (10, 6))
ax = fig.add_subplot(1, 1, 1)
ax.plot(np.real(poles), np.imag(poles), '--')
ax.plot(np.real(polesE), np.imag(polesE), 'ko', markersize = 4)
ax.set_xlabel('Real')
ax.set_ylabel('Imag')
ax.grid()
plt.show()
else:
poles = approx.getPoles()
print("Poles:\n{}".format(poles))
print("\n")

Event Timeline

c4science · Help