Commit 74161976 authored by Lars Stietz's avatar Lars Stietz
Browse files

cleaning up code and creating plots

parent 45ab9b6c
......@@ -4,14 +4,23 @@ from copy import deepcopy
def woSubstitute(u_0, Q, F):
sz = Q[1:, 1:].shape[1]
"""
Parameters:
u_0 (numpy.ndarray) : vector of initial values
Q (numpy.ndarray) : matrix with quadrature weights (node to node)
F (numpy.ndarray) : matrix assembling the source term of ODE
Returns:
X (num
"""
I = np.eye(len(u_0))
sol = la.solve(I-Q@F, u_0)
sol = la.solve(I - Q @ F, u_0)
X, V = np.split(sol, 2)
return X, V
def wSubstitute(u_0, Q, k):
def wSubstitute(u_0, Q, k, f):
sz = Q[1:, 1:].shape[1]
x_0, v_0 = np.split(u_0, 2)
X = la.solve(np.eye(sz) + k * np.linalg.matrix_power(Q[1:, 1:], 2), x_0 + Q[1:, 1:] @ v_0)
......@@ -43,48 +52,42 @@ def leapFrogSweep(x_0, v_0, dt, f):
return X[1:], V[1:]
def secondOrderSDCsweep(x,v,dt,S,f):
def secondOrderSDCsweep(x, v, dt, S, f):
x_old = deepcopy(x)
v_old = deepcopy(v)
for m in range(1,len(dt)):
x[m] = x[m-1] + dt[m]*(v[m-1]-v_old[m-1]) + 0.5*(dt[m]**2)*(f(x[m-1])-f(x_old[m-1])) + S[m+1,1:]@v_old
v[m] = v[m-1] + 0.5*dt[m]*(f(x[m-1])+f(x[m])) - 0.5*dt[m]*((f(x_old[m-1])+f(x_old[m]))) + S[m+1,1:]@f(x_old)
for m in range(1, len(dt)):
x[m] = x[m - 1] + dt[m] * (v[m - 1] - v_old[m - 1]) + 0.5 * (dt[m] ** 2) * (f(x[m - 1]) - f(x_old[m - 1])) \
+ S[m + 1, 1:] @ v_old
v[m] = v[m - 1] + 0.5 * dt[m] * (f(x[m - 1]) + f(x[m])) - 0.5 * dt[m] * ((f(x_old[m - 1]) + f(x_old[m]))) \
+ S[m + 1, 1:] @ f(x_old)
return x,v
return x, v
def secondOrderSDCsweep_opt(x, v, U_0, Qdelta, Q, k):
def secondOrderSDCsweep_opt(x, v, Qdelta, Q, k):
x_old = deepcopy(x)
v_old = deepcopy(v)
sz = len(x_old)
for m in range(1, sz):
x[m] = x[0] + np.dot(Qdelta[m, :m], (v[:m]-v_old[:m])) - k*np.dot(Qdelta[m, sz:sz+m], (x[:m]-x_old[:m])) + np.dot(Q[m, :], v_old)
v[m] = v[0] - k*np.dot(Qdelta[m+1, sz:], x) + k*np.dot(Qdelta[m, sz:], x_old)-k*np.dot(Q[m, :], x_old)
x[m] = x[0] + np.dot(Qdelta[m, :m], (v[:m] - v_old[:m])) - k * np.dot(Qdelta[m, sz:sz + m], (x[:m] - x_old[:m])) \
+ np.dot(Q[m, :], v_old)
v[m] = v[0] - k * np.dot(Qdelta[m + 1, sz:], x) + k * np.dot(Qdelta[m, sz:], x_old) - k * np.dot(Q[m, :], x_old)
return x, v
def asmFmat(k,num_nodes):
F = np.block([[np.zeros([num_nodes,num_nodes]),np.eye(num_nodes)],[-k*np.eye(num_nodes),np.zeros([num_nodes,num_nodes])]])
def asmFmat(k, num_nodes):
F = np.block([[np.zeros([num_nodes, num_nodes]), np.eye(num_nodes)],
[-k * np.eye(num_nodes), np.zeros([num_nodes, num_nodes])]])
return F
def FU(U):
fu = F@U
return fu
def secondOrderMat(X, V, U_0, Qdelta, Q, k):
x_old = deepcopy(X)
U = np.concatenate([X, V])
I = np.eye(len(U))
F = asmFmat(k, int(0.5*len(U)))
Uout = np.linalg.solve(I-Qdelta@F, U_0) + U - np.linalg.solve(I-Qdelta@F, (I-Q@F)@U)
X,V = np.split(Uout, 2)
return X,V
F = asmFmat(k, int(0.5 * len(U)))
Uout = np.linalg.solve(I - Qdelta @ F, U_0) + U - np.linalg.solve(I - Qdelta @ F, (I - Q @ F) @ U)
X, V = np.split(Uout, 2)
return X, V
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
#import matplotlib.gridspec as gridspec
# import matplotlib.gridspec as gridspec
from gauss_lobatto import CollGaussLobatto
import SDCsecondOrder as SDC
from MatrixSDC import MatrixSDC as mat
......@@ -10,40 +10,50 @@ import optimise_qdelta as opt
'''Initialize Problem'''
##################################################################################################################
k = 1
k = 1 # ODE parameter
omega = np.sqrt(k)
x_0 = 1
v_0 = 0
t_start = 0.0
t_end = 0.5
x_0 = 1 # start position
v_0 = 0 # start velocity
t_start = 0.0 # starting time
t_end = 0.5 # end time
# dx/dt = v
# dv/dt = f(x)
def f(x):
return -k*x
return -k * x # source term of ODE
##################################################################################################################
'''Analytical solution'''
##################################################################################################################
def x(t):
return x_0*np.cos(omega*t)+v_0/omega*np.sin(omega*t)
return x_0 * np.cos(omega * t) + v_0 / omega * np.sin(omega * t)
def v(t):
return -omega*x_0*np.sin(omega*t)+v_0*np.cos(omega*t)
return -omega * x_0 * np.sin(omega * t) + v_0 * np.cos(omega * t)
##################################################################################################################
'''SDC'''
##################################################################################################################
num_tsteps = 2
num_nodes = 3 #optim only tested for three nodes
num_nodes = 3 # optim only tested for three nodes
iter = 10
aim = "norm"
plot = 1
x_error = np.empty(0)
v_error = np.empty(0)
#delta_t = np.empty(0)
# delta_t = np.empty(0)
subinterval = np.linspace(t_start, t_end, num_tsteps)
#delta_t = (t_end - t_start) / num_tsteps
# delta_t = (t_end - t_start) / num_tsteps
it_steps = np.zeros(num_tsteps, dtype=int)
sol_x = np.empty(0)
......@@ -79,9 +89,9 @@ for i in range(num_tsteps - 1):
res = 1
it = 0
while res > 1e-15 and it < 15:
#X, V = SDC.secondOrderSDCsweep(Xold, Vold, coll.delta_m, coll.Smat, f)
X, V = SDC.secondOrderSDCsweep_opt(Xold, Vold, Q_delta, coll.Qmat[1:,1:], k)
#X, V = SDC.secondOrderMat(Xold, Vold, initial_val, matrices.Qdelta, matrices.Qmat, k)
# X, V = SDC.secondOrderSDCsweep(Xold, Vold, coll.delta_m, coll.Smat, f)
X, V = SDC.secondOrderSDCsweep_opt(Xold, Vold, Q_delta, coll.Qmat[1:, 1:], k)
# X, V = SDC.secondOrderMat(Xold, Vold, initial_val, matrices.Qdelta, matrices.Qmat, k)
Xopt, Vopt = SDC.secondOrderSDCsweep_opt(Xoldopt, Voldopt, opt_Qdelta, coll.Qmat[1:, 1:], k)
Xold = deepcopy(X)
Vold = deepcopy(V)
......@@ -109,7 +119,7 @@ for i in range(num_tsteps - 1):
v_0 = sol_v[-1]
kin_Energy = 0.5 * np.square(sol_v) + 0.5 * omega ** 2 * np.square(sol_x)
x_error = np.append(x_error, la.norm(x(nodes) - sol_x) / la.norm(sol_x))
v_error = np.append(v_error,la.norm(v(nodes)-sol_v)/la.norm(sol_v))
v_error = np.append(v_error, la.norm(v(nodes) - sol_v) / la.norm(sol_v))
if plot:
px = 1 / plt.rcParams['figure.dpi']
......
......@@ -18,8 +18,9 @@ matrices = mat(coll.delta_m, coll.Qmat[1:, 1:])
F = SDC.asmFmat(1, num_nodes)
Q = matrices.Qmat
Q_delta = matrices.Qdelta
aim = "norm"
opt_Qdelta,norm_val_opt = opt.get_optimal_qdelta(Q, Q_delta, F, "norm")
opt_Qdelta, norm_val_opt = opt.get_optimal_qdelta(Q, Q_delta, F, "norm")
norm_k = np.zeros(len(k))
norm_mat = np.zeros(len(k))
......@@ -37,5 +38,6 @@ ax3.set_xlabel('k')
ax3.set_ylabel('$ \|.\|$')
fig2.legend()
ax3.set_title('Iteration Matrix Properties')
plt.savefig("aim_{}.pdf".format(aim))
plt.show()
......@@ -34,12 +34,9 @@ def v(t):
##################################################################################################################
num_tsteps = 10
num_nodes = 3 #optim only tested for three nodes
iter = 10
aim = "norm"
plot = 1
x_error = np.empty(0)
v_error = np.empty(0)
#delta_t = np.empty(0)
subinterval = np.linspace(t_start, t_end, num_tsteps)
......@@ -72,8 +69,8 @@ spec_opt_Qdelta, spec_opt_val = opt.get_optimal_qdelta(Q, Q_delta, F, "spec")
spec_opt_norm = la.norm(I-la.inv(I-spec_opt_Qdelta@F)@(I-Q@F),np.inf)
norm_opt_Qdelta, norm_opt_val = opt.get_optimal_qdelta(Q, Q_delta, F, "norm")
norm_opt_spec = max(abs(la.eigvals(I-la.inv(I-norm_opt_Qdelta@F)@(I-Q@F))))
rowname = [r"$\rho(E)$", r"$\| E \|$"]
colname = ["Without optimisation", r"Optimised wrt $\rho(E)$", r"Optimised wrt $\| E \|$"]
rowname = [r"$\rho(\mathbf{E})$", r"$\Vert \mathbf{E} \Vert$"]
colname = ["Without optimization", r"Optimized with $C_\mathrm{spec}$", r"Optimized with $C_\mathrm{norm}$"]
cell_data =[[spec_val, spec_opt_val, norm_opt_spec], [norm_val, spec_opt_norm, norm_opt_val]]
X_0 = np.ones(num_nodes) * x_0
......@@ -118,7 +115,6 @@ for i in range(num_tsteps - 1):
res_norm_norm = np.append(res_norm_norm, la.norm(res_x_norm, np.inf))
res = la.norm(res_x_spec, np.inf)
it = it + 1
print(it)
it_steps[i + 1] = it
sol_x = np.append(sol_x, X)
......@@ -129,8 +125,8 @@ for i in range(num_tsteps - 1):
sol_v_norm = np.append(sol_v_norm, Vnorm)
res_x = X_0 + coll.Qmat[1:, 1:] @ V - X
res_x_spec = X_0 + coll.Qmat[1:, 1:] @ Vspec - Xspec
res_x_norm = X_0 + coll.Qmat[1:, 1:] @ Vnorm - Xnorm
res_x_spec = X_0_spec + coll.Qmat[1:, 1:] @ Vspec - Xspec
res_x_norm = X_0_norm + coll.Qmat[1:, 1:] @ Vnorm - Xnorm
res_it_x = np.append(res_it_x, res_x)
res_it_x_spec = np.append(res_it_x_spec, res_x_spec)
res_it_x_norm = np.append(res_it_x_norm, res_x_norm)
......@@ -173,50 +169,60 @@ if plot:
fax = fig.add_subplot(gs[2, 4:])
ax1.plot(sol_x_spec, sol_v_spec)
ax1.set_title(r"SDC solution with optimisation wrt $\rho(E)$", usetex = True)
#ax1.set_title(r"SDC solution with optimisation wrt $\rho(E)$", usetex = True)
ax1.set_ylabel("v(t)")
ax1.set_xlabel("x(t)")
#plt.savefig("solution_k-{}.pdf".format(k))
ax2.semilogy(nodes, abs(res_it_x), label="Res of x(t)")
ax2.semilogy(nodes, abs(res_it_x_spec), label=r"Res of x(t) optimised wrt $\rho(E)$")
ax2.semilogy(nodes, abs(res_it_x_norm), label=r"Res of x(t) optimised wrt $\|E\|$")
ax2.set_title("Residual")
ax2.semilogy(nodes, abs(res_it_x_spec), label=r"Res of x(t) optimized wrt $\rho(E)$")
ax2.semilogy(nodes, abs(res_it_x_norm), label=r"Res of x(t) optimized wrt $\Vert E \Vert$")
#ax2.set_title("Residual")
ax2.set_xlabel("t")
ax2.legend(loc="upper center")
fax.semilogy(res_norm_x, label="Standard Residual")
fax.semilogy(res_norm_spec, label=r"wrt $\rho(E)$ optimised")
fax.semilogy(res_norm_norm, label=r"wrt $\|E\|$ optimised")
fax.semilogy(res_norm_spec, label=r"wrt $\rho(E)$ optimized")
fax.semilogy(res_norm_norm, label=r"wrt $\Vert E \Vert$ optimized")
# fax.set_ylim(0.49, 0.51)
fax.set_ylabel(r"$|\mathrm{res}|$", usetex=True)
fax.set_xlabel("#iter")
fax.set_title("Residual norm after each iteration")
#fax.set_title("Residual norm after each iteration")
fax.legend()
#plt.savefig('k-{}_tsteps-{}.pdf'.format(k, num_tsteps))
plt.figure()
plt.semilogy(res_norm_x[it_steps[0]:(it_steps[1] + it_steps[2])], label="Standard Residual")
plt.semilogy(res_norm_spec[it_steps[0]:(it_steps[1] + it_steps[2])], label=r"Optimised Residual wrt $\rho(E)$")
plt.semilogy(res_norm_norm[it_steps[0]:(it_steps[1] + it_steps[2])], label=r"Optimised Residual wrt $\| \mathrm{E} \|$")
#plt.semilogy(np.power(spec_opt_val, np.tile(range(0,it_steps[1]),2)), label = "spec(E)^k")
#plt.semilogy(np.power(norm_opt_val, np.tile(range(0, it_steps[1]), 2)), label="norm(E)^k")
# fax. set_ylim(0.49, 0.51)
plt.ylabel(r"$| \mathrm{res} |$", usetex=True)
#plt.xlabel("#iter")
plt.title("Residual norm after each iteration for two timesteps")
plt.semilogy(res_norm_x[it_steps[0]:(it_steps[1] + it_steps[2])], label="no optimization")
plt.semilogy(res_norm_spec[it_steps[0]:(it_steps[1] + it_steps[2])], label=r" $C_\mathrm{spec}$")
plt.semilogy(res_norm_norm[it_steps[0]:(it_steps[1] + it_steps[2])], label=r"$C_\mathrm{norm}$")
plt.ylabel(r"$\Vert \textbf{R}_{x} \Vert$", usetex=True)
plt.xlabel("#iterations")
#plt.title("Residual norm after each iteration for two timesteps")
plt.legend()
plt.xticks([])
the_table = plt.table(cellText=cell_data,
rowLabels=rowname,
colLabels=colname,
loc='bottom')
loc='top')
plt.tight_layout()
plt.savefig('res_k-{}_tsteps-{}.pdf'.format(k, num_tsteps))
plt.figure()
plt.plot(nodes,kin_Energy, label = "standard")
plt.plot(nodes, kin_Energy_spec, '--', label = "spec opt")
plt.plot(nodes, kin_Energy_norm, label = "norm opt")
plt.ylabel("total energy")
plt.plot(nodes, kin_Energy, label="no optimization")
plt.plot(nodes, kin_Energy_spec, '--', label =r" $C_\mathrm{spec}$")
plt.plot(nodes, kin_Energy_norm, label =r"$C_\mathrm{norm}$")
plt.ylabel(r"$E(t)$")
plt.xlabel("time [t]")
plt.title("Total energy calculated by different SDC optimized iterations")
#plt.title("Total energy calculated by different SDC optimized iterations")
plt.legend()
plt.savefig('energy_k-{}_tsteps-{}.pdf'.format(k, num_tsteps))
plt.show()
plt.axes().set_aspect('equal')
plt.plot(sol_x_spec, sol_v_spec)
# ax1.set_title(r"SDC solution with optimisation wrt $\rho(E)$", usetex = True)
plt.ylabel("v(t)")
plt.xlabel("x(t)")
#plt.savefig("solution_k-{}.pdf".format(k))
\ No newline at end of file
import numpy as np
from numpy import linalg as LA
import scipy.optimize as opt
from pseudo_spectral_radius import pseudo_spectral_radius
import scipy.sparse as sparse
"""Optimization with 3 nodes explicitly"""
def objective_func(param, Q, F, aim):
def objective_func(param, Q, F, aim, eps = 1e-5):
"""
Objective function to be minimized.
......@@ -24,6 +25,10 @@ def objective_func(param, Q, F, aim):
return max(abs(LA.eigvals(I-LA.inv(I-Q_delta@F)@(I-Q@F))))
elif aim == "norm":
return LA.norm(I-LA.inv(I-Q_delta@F)@(I-Q@F), np.inf)
elif aim == "pseudo":
E = pseudo_spectral_radius(Q_delta, eps)
sol, _, _, _ = E.get_psr()
return sol
def constraint_func(param, Q, F):
......@@ -62,8 +67,8 @@ def get_optimal_qdelta(Q, Q_delta, F, aim):
param = assemble_param(Q_delta)
f = lambda a: objective_func(a, Q, F, aim)
con = lambda a: constraint_func(a, Q, F)
nlc = opt.NonlinearConstraint(con, [0, 0], [0.5, 1])
# con = lambda a: constraint_func(a, Q, F)
# nlc = opt.NonlinearConstraint(con, [0, 0], [0.5, 1])
print("Objective function before optimisation: %5.3e" % objective_func(param, Q, F, aim))
print("\n ... optimizing ... \n")
# sol = opt.minimize(f, param, constraints = nlc, method='SLSQP', options={'maxiter': 10000, 'fatol':1e-16})
......
from scipy.optimize import NonlinearConstraint
from scipy.optimize import minimize
import scipy.sparse as sparse
from scipy.linalg import svdvals
import numpy as np
'''
Computes the epsilon pseudo spectral radius of the matrix E.
'''
class pseudo_spectral_radius(object):
def __init__(self, E, eps):
assert isinstance(eps, float) and eps > 0, "Parameter eps must be a positive real number"
self.eps = eps
self.E = E
self.n = np.shape(E)[0]
assert self.n==np.shape(E)[1], "Matrix E must be square"
self.Id = np.eye(self.n)
'''
The constraint makes sure we are searching along an epsilon isoline in the pseudo-spectrum
'''
def constraint(self, x):
z = x[0] + 1j*x[1]
# See algorithm on p. 371, Chapter 39 in the Trefethen book
M = z*self.Id - self.E
sv = svdvals(M)
return np.min(sv)
'''
Minimise 1/||x||_2 while the constraint keeps us on the isoline with maximise distance from the origin and yield the spectral radius.
'''
def target(self, x):
return 1.0/np.linalg.norm(x, 2)**2
'''
Computes the epsilon-pseudospectral radius of the matrix E.
'''
def get_psr(self):
# The constraint will keep us on an isoline where the minimum singular value of z I - E equals eps, meaning that || (z I - E)^(-1) ||_2 = 1\eps - we allow for a bit of variation (1e-9) to help with convergence of the optimiser
nlc = NonlinearConstraint(self.constraint, self.eps-1e-9, self.eps+1e-9)
# Now run the minimiser to minimise 1/||x||_2^2 while keeping min(svd(zI-E))=eps
result = minimize(self.target, [np.sqrt(self.eps), np.sqrt(self.eps)], constraints=nlc, tol=1e-10, method='trust-constr', options={'xtol': 1e-10, 'gtol': 1e-10, 'maxiter': 500})
# The eps - pseudo spectral radius corresponds to the maximum distance from the origin on the eps isoline of min(svd(zI-E))
return np.linalg.norm(result.x,2), result.x, self.target(result.x), self.constraint(result.x)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment