Hello,
I am a bit confused about the implementation of SDIRK method in ngsolve when i studied the SDIRK method in i-tutorial unit 3.1
introduction of SDIRK method
Let us consider more sophisticated time integration methods, SDIRK methods. To simplify presentation we set f=0.
SDIRK methods for the semi-discrete problem
\frac{d}{dt} u = M^{-1} F(u) = -M^{-1} \cdot A u
are of the form:
u^{n+1} = u^n + \Delta t M^{-1} \sum_{i=0}^{s-1} b_i k_i
with k_i = - A u_i , where u_i is the solution to
(M + a_{ii} \Delta t A) u_i = M u^n - \Delta t \sum_{i=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.
The coefficients a,b and c are stored in the butcher tableau:
\begin{array}{c|cccc}
c_0 & a_{00} & 0 & \ddots& \\
c_1 & a_{10} & a_{11} & 0 & \ddots \\
\vdots & \vdots & \vdots & \ddots & 0\\
c_{s-1} & a_{s-1,0} & a_{s-1,1} & \dots& a_{s-1,s-1} \\ \hline
& b_{0} & b_1 & \dots& b_{s-1} \\
\end{array}
For an SDIRK method we have a^\ast = a_{ii},~i=0,..,s-1.
implementation
- butcher tableau
class sdirk5: #order 4
stages = 5
a = [[1/4, 0, 0, 0, 0],
[1/2, 1/4, 0, 0, 0],
[17/50,-1/25, 1/4, 0, 0],
[371/1360, -137/2720, 15/544, 1/4,0],
[25/24, -49/48, 125/16, -85/12, 1/4]]
b = [25/24, -49/48, 125/16, -85/12, 1/4]
c = [1/4, 3/4, 11/20, 1/2, 1]
astar = 1/4
- time stepping function
def TimeStepping_app4(m, a, butchertab, dt,
initial_cond = None, t0 = 0, tend = 2,
nsamples = 10):
mstar.AsVector().data = m.mat.AsVector() + butchertab.astar * dt * a.mat.AsVector()
invmstar = mstar.Inverse(freedofs=fes.FreeDofs())
invmass = m.mat.Inverse(freedofs=fes.FreeDofs())
rhsi, Mu0, ui = gfu.vec.CreateVector(),gfu.vec.CreateVector(),gfu.vec.CreateVector()
k = [gfu.vec.CreateVector() for i in range(butchertab.stages)]
if initial_cond:
gfu.Set(initial_cond)
cnt = 0; time = t0
sample_int = int(floor(tend / dt / nsamples)+1)
gfuD = GridFunction(gfu.space)
gfut = GridFunction(gfu.space,multidim=0)
gfut.AddMultiDimComponent(gfu.vec)
while time < tend - 0.5 * dt:
Mu0.data = m.mat * gfu.vec
for i in range(butchertab.stages):
rhsi.data = Mu0
for j in range(0,i):
rhsi.data += dt * butchertab.a[i][j] * k[j]
t.Set(time+butchertab.c[i]*dt)
gfu.Set(uD,BND)
ui.data = gfu.vec; rhsi.data -= mstar * ui
ui.data += invmstar * rhsi
k[i].data = - a.mat * ui
t.Set(time+dt)
gfu.Set(uD,BND)
Mu0.data -= m.mat * gfu.vec
for i in range(0,butchertab.stages):
Mu0.data += dt * butchertab.b[i] * k[i]
gfu.vec.data += invmass * Mu0
print("\r",time,end="")
if cnt % sample_int == 0:
gfut.AddMultiDimComponent(gfu.vec)
cnt +=
questions
- inconsistency between formulas and codes
when we use the results of previous sub-steps u_{j} to compute u_{i}, i think the formula should be of the following form (formula in the tutotial made a small mistake in \Sigma)
(M + a_{ii} \Delta t A) u_i = M u^n - \Delta t \sum_{j=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.
The implementation of the formula is as follows (i add some comments with #====)
Mu0.data = m.mat * gfu.vec
for i in range(butchertab.stages):
#====now we have rhs = Mu^{n}====
rhsi.data = Mu0
for j in range(0,i):
#====why not rhsi.data -= dt * butchertab.a[i][j] * k[j]====
rhsi.data += dt * butchertab.a[i][j] * k[j]
t.Set(time+butchertab.c[i]*dt)
#====time-dependent boundary conditions====
gfu.Set(uD,BND)
ui.data = gfu.vec
rhsi.data -= mstar * ui
#====homogeneous + BC====
ui.data += invmstar * rhsi
k[i].data = - a.mat * ui
t.Set(time+dt)
Is there a sign problem? (rhsi.data += dt * butchertab.a[i][j] * k[j] for M u^n - \Delta t \sum_\limits{j=0}^{i-1} a_{ij} k_j, \quad i=0,..,s-1.)