Commit 45ab9b6c authored by Lars Stietz's avatar Lars Stietz
Browse files

added comments

parent 18590e74
import numpy as np
class MatrixSDC:
def __init__(self,dt,Q):
"""
A class to assemble matrices for the SDC method for second order ODEs
The Matrices created here are according to the paper DOI: 10.1016/j.jcpx.2019.100036
...
Attributes:
dt (numpy.ndarray): array of node to node distances in paper delta tau
Q (numpy.ndarray): matrix containing the quadrature weights for zero to node 1D case
Qmat (numpy.ndarray): matrix containing the quadrature weights for zero to node 2D case
Qe (numpy.ndarray): matrix in paper Q_{delta,E}
Qi (numpy.ndarray): matrix in paper Q_{delta,I}
Qt (numpy.ndarray): matrix in paper Q_{delta,T}
Qdelta (numpy.ndarray): matrix in paper Q_{delta}
"""
def __init__(self, dt, Q):
"""
Parameters
----------
dt : numpy.ndarray
Array of distances between quadrature nodes
Q : numpy.ndarray
Matrix with quadrature weights (zero to node)
"""
#set node to node distances and quadrature weights
self.dt = dt
self.Q = Q
#
self.Qmat = self._getQmat
self.Qe = self._getQe
self.Qi = self._getQi
self.Qq = self._getQq
self.Qt = self._getQt
self.Qdelta = self._getQdelta
@property
def _getQmat(self):
Qmat = np.kron(np.eye(2),self.Q)
"""
Returns:
numpy.ndarray: matrix containing the quadrature weights for 2D case (position, velocities)
"""
Qmat = np.kron(np.eye(2), self.Q)
return Qmat
@property
def _getQe(self):
"""
Returns:
numpy.ndarray: matrix containing the node to node distances (dt2,...dtM,0)
in a strictly lower triangular matrix
"""
sz = len(self.dt)
Qe = np.zeros([sz, sz])
for i in range(1, sz):
......@@ -25,29 +61,31 @@ class MatrixSDC:
@property
def _getQi(self):
"""
Returns:
numpy.ndarray: matrix containing the node to node distances (dt1,...,dtM) in a lower triangular matrix
"""
sz = len(self.dt)
Qi = np.zeros([sz,sz])
Qi = np.zeros([sz, sz])
for i in range(sz):
Qi[i,:i+1] = self.dt[:i+1]
Qi[i, :i+1] = self.dt[:i+1]
return Qi
@property
def _getQq(self):
Qq = 0.5*np.square(self.Qe)
return Qq
@property
def _getQt(self):
"""
Returns:
numpy.ndarray: matrix assembled form 0.5(Qe+Qi)
"""
Qt = 0.5*(self.Qe+self.Qi)
return Qt
@property
def _getQdelta(self):
"""
Returns:
numpy.ndarray: matrix needed for SDC iteration
"""
sz = len(self.dt)
Qdelta = np.block([[self.Qe,self.Qq],[np.zeros([sz,sz]),self.Qt]])
Qdelta = np.block([[self.Qe, 0.5*np.square(self.Qe)], [np.zeros([sz, sz]), self.Qt]])
return Qdelta
......@@ -25,7 +25,7 @@ dt = np.diff(subinterval)
# initialize Vector
u_0 = u_0M = 1
l = -5j
l = -10j
"""
......@@ -66,7 +66,7 @@ if solving == 1:
res_2 = 2
res_2M = 2
while res_2 > 1e-15:
u = sweep(u,coll.delta_m,coll.Smat,coll.Qmat,f)
u = sweepMat(u_Mat,coll.Qmat,asmFastMat(coll.delta_m),f)
r = res(u,coll.Qmat,u_0,f)
res_2 = la.norm(r)
res_it = np.append(res_it,res_2)
......@@ -80,9 +80,9 @@ if solving == 1:
# standard Q_fast
#Q_fast = asmFastMat(coll.delta_m,dt[i])
u_Mat = sweepMat(u_Mat,coll.delta_m,dt[i],coll.Qmat,Q_fast,f)
u_Mat = sweepMat(u_Mat,coll.Qmat,Q_fast,f)
r_Mat = res(u_Mat,coll.Qmat,u_0,f)
r_Mat = res(u_Mat,coll.Qmat,u_0M,f)
res_2M = la.norm(r_Mat)
res_itM = np.append(res_itM,res_2M)
......
......@@ -3,14 +3,14 @@ import scipy.linalg as la
from copy import deepcopy
def woSubstitute(u_0, Q, k):
def woSubstitute(u_0, Q, F):
sz = Q[1:, 1:].shape[1]
collMat = np.block([[np.eye(sz), -Q[1:, 1:]], [k * Q[1:, 1:], np.eye(sz)]])
sol = la.solve(collMat, u_0)
I = np.eye(len(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):
sz = Q[1:, 1:].shape[1]
x_0, v_0 = np.split(u_0, 2)
......@@ -18,6 +18,7 @@ def wSubstitute(u_0, Q, k):
V = v_0 + Q[1:, 1:] @ f(X)
return X, V
def verletSweep(x_0, v_0, dt, f):
X = np.zeros(dt.shape[0] + 1)
X[0] = x_0
......@@ -28,6 +29,7 @@ def verletSweep(x_0, v_0, dt, f):
V[i + 1] = V[i] + 0.5 * dt[i] * (f(X[i]) + f(X[i + 1]))
return X[1:], V[1:]
def leapFrogSweep(x_0, v_0, dt, f):
X = np.zeros(dt.shape[0] + 1)
X[0] = x_0
......@@ -40,26 +42,46 @@ def leapFrogSweep(x_0, v_0, dt, f):
V[i + 1] = vhalf + 0.5 * dt[i] * f(X[i + 1])
return X[1:], V[1:]
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)
print(max(abs(x_old-x)))
return x,v
def secondOrderSDCsweep_opt(x, v, U_0, 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)
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])]])
return F
def FU(U):
fu = F@U
return fu
def secondOrderMat(X,V,U_0,Qdelta,Q,k):
U = np.concatenate([X,V])
F = asmFmat(k,int(0.5*len(U)))
Uout = np.linalg.solve(np.eye(len(U))-Qdelta@F,U_0)+U - np.linalg.solve(np.eye(len(U))-Qdelta@F,(np.eye(len(U))-Q@F)@U)
X,V = np.split(Uout,2)
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
......
import numpy as np
from copy import deepcopy
"""
Assemble all matrices needed for SDC for first order ODEs according to DOI: 10.1137/16M1060078
# matrix assembling routine
def asmSlowMat(dt):
sz = len(dt)
Q = np.zeros([sz,sz])
for i in range(1,sz):
Q[i,:i] = dt[1:i+1]
The functions assemble the matrices needed for SDC with routines for sweep or matrix formulation of one SDC iteration.
Only considered quadrature nodes, where first node is same as tstart.
"""
# adjustment for first node not equal to tstart needed
def asmSlowMat(d_tau):
"""
Assembles Q_delta matrix for slow part from vector dt
Parameters:
d_tau (numpy.ndarray) : vector containing the node to node distances
Returns:
Q (numpy.ndarray) : strictly lower triangular matrix corresponding to Q^{slow}_{\delta}
"""
sz = len(d_tau)
Q = np.zeros([sz, sz])
for i in range(1, sz):
Q[i, :i] = d_tau[1:i+1]
return Q
def asmFastMat(dt):
sz = len(dt)
Q = np.zeros([sz,sz])
for i in range(0,sz):
Q[i,:i+1] = dt[:i+1]
def asmFastMat(d_tau):
"""
Assembles Q_delta matrix for fast part from vector d_tau
Parameters:
d_tau (numpy.ndarray) : vector containing the node to node distances
Returns:
Q (numpy.ndarray) : lower triangular matrix corresponding to Q^{fast}_{\delta}
"""
sz = len(d_tau)
Q = np.zeros([sz, sz])
for i in range(0, sz):
Q[i, :i+1] = d_tau[:i+1]
return Q
# function for SDC sweep in one iteration step
def sweep(u,d_tau,Smat,Qmat,f):
def sweep(u, d_tau, Smat, f):
"""
adjustment for strt/nd=False still needed
Sweep formutlation of one SDC iteration step.
Parameters:
u (numpy.ndarray) : vector of ODE solution before one SDC iteration step (sweep)
d_tau (numpy.ndarray) : vector containing the node to node distances
Smat (numpy.ndarray) : matrix containing the quadrature weights for node to node
f (function) : source term of the ODE
Returns:
u (numpy.ndarray) : vector ODE solution after one SDC iteration step (sweep)
"""
num_nodes = len(u)-1
u_old = deepcopy(u)
quad = f(1)*Smat[1:,1:]@u_old[1:-1]
for m in range(1,num_nodes):
quad = f(1)*Smat[1:, 1:]@u_old[1:-1]
for m in range(1, num_nodes):
u[m] = 1/(1-f(1)*d_tau[m-1])*(u[m-1]-d_tau[m-1]*f(u_old[m]) + quad[m-1])
return u
# function for SDC sweep with matrices
def sweepMat(u,d_tau,dt,Qmat,Qs,f):
def sweepMat(u, Qmat, Qs, f):
"""
Matrix formulation of one SDC iteration step
Parameters:
u (numpy.ndarray) : vector of ODE solution before one SDC iteration step (matrix formulation)
Qmat (numpy.ndarray) : matrix containing the quadrature weights for zero to node
Qs (numpy.ndarray) : matrix corresponding to Q^{fast}_{\delta}
f (function) : source term of the ODE
Returns:
u (numpy.ndarray) : vector of ODE solution after one SDC iteration step (matrix formulation)
-------
"""
u_old = deepcopy(u)
RHS = u[0]-dt*f(1)*Qs@u_old[1:-1]+f(1)*Qmat[1:,1:]@u_old[1:-1]
RHS = u[0]-f(1)*Qs@u_old[1:-1]+f(1)*Qmat[1:, 1:]@u_old[1:-1]
LHS = np.eye(np.shape(Qs)[0])-f(1)*Qs
u[1:-1] = np.linalg.solve(LHS,RHS)
u[1:-1] = np.linalg.solve(LHS, RHS)
return u
# function to calculate the residual
def res(u, Qmat,u_0,f):
def res(u, Qmat, u_0, f):
"""
adjustment for strt/nd=False still needed
Calculates the residual
Parameters:
u (numpy.ndarray) : vector of the solution of the ODE
Qmat (numpy.ndarray) : matrix containing the quadrature weights for zero to node
u_0 (numpy.ndarray) : vector of initial value
f (function) : source term of the ODE
Returns:
r (numpy.ndarray) : residual of the solution
"""
r = u_0 + Qmat[1:,1:]@f(u[1:-1])-u[1:-1]
return r
\ No newline at end of file
r = u_0 + Qmat[1:, 1:]@f(u[1:-1])-u[1:-1]
return r
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