Commit de837464 authored by Davide Torlo's avatar Davide Torlo
Browse files

added a reference to the original RRK work in the notebook

parent 452aa909
%% Cell type:markdown id: tags:
 
# Relaxation Deferred Correction notebook
%% Cell type:markdown id: tags:
This notebook is an extension of the original [Relaxation Runge--Kutta Notebook](https://github.com/ketch/RRK_rr) provided by D. Ketcheson and H. Ranocha on the work of H. Ranocha, M. Sayyari, L. Dalcin, M. Parsani and D. I. Ketcheson on [Relaxation Runge--Kutta Methods: Fully Discrete Explicit Entropy-Stable Schemes for the Compressible Euler and Navier--Stokes Equations](https://doi.org/10.1137/19M1263480).
Here we use the Deferred Correction method as a Runge--Kutta method. The implementation of the Deferred Correction method as a Runge--Kutta method is included in [DeC.py](DeC.py) and explained in details in the work [Relaxation Deferred Correction Methods and their Applications to Residual Distribution Schemes ](https://arxiv.org/search/?searchtype=author&query=Torlo%2C+D) by R. Abgrall, E. Le Mélédo, P. Öffner and D. Torlo.
%% Cell type:markdown id: tags:
# Setup
 
Note that this notebook was developed with NodePy version 0.7.
 
%% Cell type:code id: tags:
 
``` python
import numpy as np
import matplotlib.pyplot as plt
from nodepy import rk, stability_function
from DeC import *
 
rk4 = rk.loadRKM('RK44').__num__()
rk4x2 = rk4*rk4
ssp2 = rk.loadRKM('SSP22').__num__()
ssp3 = rk.loadRKM('SSP33').__num__()
ssp104 = rk.loadRKM('SSP104').__num__()
merson4 = rk.loadRKM('Merson43').__num__()
bs5 = rk.loadRKM('BS5').__num__()
 
decEqRKs =[]
decGLBRKs =[]
for order in range(2,10):
A,b,c=compute_RK_from_DeC(order-1,order,"equispaced") # subtimesteps, iterations, nodetype{"equispaced","gaussLobatto"}
decEqRKs.append(rk.ExplicitRungeKuttaMethod(A,b))
decEqRKs[order-2].name="DeC%dEq"%(order)
A,b,c=compute_RK_from_DeC(np.int(np.ceil(order/2.)),order,"gaussLobatto") # subtimesteps, iterations, nodetype{"equispaced","gaussLobatto"}
decGLBRKs.append(rk.ExplicitRungeKuttaMethod(A,b))
decGLBRKs[order-2].name="DeC%dGLB"%(order)
 
dec3 = decEqRKs[1]
 
dec4equi = decEqRKs[2]
 
dec4galo = decGLBRKs[2]
 
 
 
#short_methods = [ssp2, ssp3, rk4, bs5, dec3]
short_methods = [ ssp2, ssp3, rk4, dec3, dec4equi, dec4galo]
#short_methods = [ dec3, dec4equi, dec4galo]
short_method_labels = ["SSPRK(2,2)=DeC(2)","SSPRK(3,3)", "RK(4,4)", "Dec3", "Dec4equi", "Dec4GaLo"]
#short_method_labels = ["Dec3", "Dec4equi", "Dec4GaLo"]
#short_method_labels = ["SSPRK(2,2)=DeC(2)", "SSPRK(3,3)", "RK(4,4)", "Dec3", "Dec4equi"]
import os
if not os.path.isdir('./figures'):
os.mkdir('./figures')
```
 
%% Cell type:code id: tags:
 
``` python
font = {'family' : 'serif',
'weight' : 'normal',
'size' : 25}
 
import matplotlib
matplotlib.rc('font', **font)
```
 
%% Cell type:code id: tags:
 
``` python
# line cycler adapted to colourblind people for only a few 10
from cycler import cycler
colors = ['#E69F00', '#56B4E9', '#009E73', '#0072B2', '#D55E00', '#CC79A7', '#F0E442']
linestyles = ['-', '--', '-.', ':', "-", "--", "-."]
markers = ['o','D','X','s','v']
colourblind_cycler = (cycler(color=colors) +
cycler(linestyle=linestyles))
plt.rc("axes", prop_cycle=colourblind_cycler)
 
plt.rc("text", usetex=True)
plt.rc("text.latex",
preamble=r"\usepackage{newpxtext}\usepackage{newpxmath}\usepackage{commath}\usepackage{mathtools}")
#plt.rc("font", family="serif", size=14.)
#plt.rc("savefig", dpi=200)
plt.rc("legend", fontsize="medium", fancybox=True, framealpha=0.5)
```
 
%% Cell type:code id: tags:
 
``` python
# line cycler adapted to colourblind people for all
from cycler import cycler
colors = ['#E69F00', '#56B4E9', '#009E73', '#0072B2', '#D55E00', '#CC79A7', '#F0E442','#F0E442']
linestyles = ['-', '--', '-.', ':', "-", "--", "-.", "--"]
markers = ['o','D','X','s','v', 'p','1','2','3']
colourblind_cycler = (cycler(color=colors) +
cycler(linestyle=linestyles))
plt.rc("axes", prop_cycle=colourblind_cycler)
 
plt.rc("text", usetex=True)
plt.rc("text.latex",
preamble=r"\usepackage{newpxtext}\usepackage{newpxmath}\usepackage{commath}\usepackage{mathtools}")
#plt.rc("font", family="serif", size=14.)
#plt.rc("savefig", dpi=200)
plt.rc("legend", fontsize="medium", fancybox=True, framealpha=0.5)
```
 
%% Cell type:code id: tags:
 
``` python
def RRK(rkm, dt, f, w0=[1.,0], t_final=1., relaxation=True,
rescale_step=True, debug=False, gammatol=0.1, print_gamma=False,
one_step=False, dissip_factor=1.0):
"""
Relaxation Runge-Kutta method implementation.
 
Options:
 
rkm: Base Runge-Kutta method, in Nodepy format
dt: time step size
f: RHS of ODE system
w0: Initial data
t_final: final solution time
relaxation: if True, use relaxation method. Otherwise, use vanilla RK method.
rescale_step: if True, new time step is t_n + \gamma dt
debug: output some additional diagnostics
gammatol: Fail if abs(1-gamma) exceeds this value
 
"""
w = np.array(w0)
t = 0
# We pre-allocate extra space because if rescale_step==True then
# we don't know exactly how many steps we will take.
ww = np.zeros([len(w0),int((t_final-t)/dt*2.5)+10000])
ww[:,0] = w.copy()
tt = [t]
ii = 0
s = len(rkm)
b = rkm.b
y = np.zeros((s,len(w0)))
max_gammam1 = 0.
gams = []
 
while t < t_final:
if t + dt >= t_final:
dt = t_final - t # Hit final time exactly
 
for i in range(s):
y[i,:] = w.copy()
for j in range(i):
y[i,:] += rkm.A[i,j]*dt*f(y[j,:])
 
F = np.array([f(y[i,:]) for i in range(s)])
 
if relaxation:
numer = 2*sum(b[i]*rkm.A[i,j]*np.dot(F[i],F[j]) \
for i in range(s) for j in range(s))
denom = sum(b[i]*b[j]*np.dot(F[i],F[j]) for i in range(s) for j in range(s))
if denom != 0:
gam = numer/denom
else:
gam = 1.
else: # Use standard RK method
gam = 1.
 
if print_gamma:
print(gam)
 
if np.abs(gam-1.) > gammatol:
print(gam)
raise Exception("The time step is probably too large.")
 
w = w + dissip_factor*gam*dt*sum([b[j]*F[j] for j in range(s)])
if (t+dt < t_final) and rescale_step:
t += dissip_factor*gam*dt
else:
t += dt
ii += 1
tt.append(t)
ww[:,ii] = w.copy()
if debug:
gm1 = np.abs(1.-gam)
max_gammam1 = max(max_gammam1,gm1)
gams.append(gam)
 
if one_step:
return w, gam
 
if debug:
print(max_gammam1)
return tt, ww[:, :ii+1], np.array(gams)
else:
return tt, ww[:,:ii+1]
```
 
%% Cell type:markdown id: tags:
 
# Introductory example (Ranocha's problem)
 
%% Cell type:code id: tags:
 
``` python
def f_prob1(u):
E = u[0]**2 + u[1]**2
return np.array([-u[1],u[0]])/E
```
 
%% Cell type:code id: tags:
 
``` python
plt.rc("font", size=20)
```
 
%% Cell type:code id: tags:
 
``` python
#Here, I have to run it again on an linux distribution! On my mac is need to much memory
u0 = np.array([1.,0.])
dt = 0.1
t_final = 10.
 
for i, method in enumerate(short_methods):
tt, uu = RRK(method, dt, f_prob1, w0=u0, t_final=t_final,
relaxation=0, debug=False)
plt.semilogy(tt[1:],uu[0,1:]**2+uu[1,1:]**2-1., lw=3, label=short_method_labels[i]);
 
 
plt.xlabel('$t$');
plt.ylabel('$\|u(t)\|^2-\|u(0)\|^2$')
plt.tight_layout()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.savefig('./figures/problem1_energy_RK.pdf')
plt.show()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
u0 = np.array([1.,0.])
dt = 0.1
 
for method in short_methods:
tt, uu = RRK(method, dt, f_prob1, w0=u0, t_final=10.,
relaxation=1, debug=False)
plt.plot(tt,uu[0,:]**2+uu[1,:]**2-1.,lw=2.5);
 
 
plt.xlabel('$t$');
plt.ylabel('$\|u(t)\|^2-\|u(0)\|^2$')
plt.tight_layout()
plt.savefig('./figures/problem1_energy_RRK.pdf')
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
plt.rc("font", size=16)
 
#names = ['SSP(2,2)', 'SSP(3,3)', 'RK(4,4)', 'BS(8,5)', 'dec_3']
dts = np.array([0.9, 0.5, 0.2, 0.1, 0.03, 0.02, 0.011, 0.005])
dts = np.logspace(-2.5,0.,11)
u0 = np.array([1.,0.])
t_final = 10.
 
fig=plt.figure(figsize=(6,6))
ax = plt.subplot(111)
lines = []
 
for j, rkm in enumerate(short_methods):
for relax in (False, True):
sols = []
for dt in dts:
 
tt, uu = RRK(rkm, dt,f_prob1, w0=u0, t_final=t_final, relaxation=relax, rescale_step=True, gammatol=1.0)
sols.append(uu[0,-1])
 
errors = []
ref = np.cos(t_final)
for i in range(len(sols)):
errors.append(abs(ref-sols[i]))
 
if relax:
plt.loglog(dts, errors,linestyle='dashed',color=colors[j],marker=markers[j],label=None);
else:
lines.append(plt.loglog(dts, errors,linestyle='-',color=colors[j],
marker=markers[j],label=short_method_labels[j])[0])
 
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
 
# Put a legend to the right of the current axis
ax.legend(lines,short_method_labels,loc='center left', bbox_to_anchor=(1, 0.5))
 
plt.xlabel('$\Delta t$')
plt.ylabel('Error');
plt.ylim([1.e-14, 10]);
 
plt.plot([1e-2,4e-2],[3.e-4, 4**2*3.e-4],'-ok',alpha=0.5)
plt.plot([1e-2,4e-2],[5.e-6, 4**3*5.e-6],'-ok',alpha=0.5)
plt.plot([1e-2,4e-2],[3.e-10, 4**4*3.e-10],'-ok',alpha=0.5)
#plt.plot([1e-2,4.e-2],[5.e-13, 4**5*5.e-13],'-ok',alpha=0.5)
#plt.plot([4e-2,16e-2],[5.e-13, 4**6*5.e-13],'-ok',alpha=0.5)
 
plt.text(1.e-2,3.e-4,'$p=2$',verticalalignment='bottom',horizontalalignment='right')
plt.text(1.e-2,5.e-6,'$p=3$',verticalalignment='top',horizontalalignment='left')
plt.text(3.e-2,6.e-9,'$p=4$',verticalalignment='bottom',horizontalalignment='left')
#plt.text(1.e-2,5.e-13,'$p=5$',verticalalignment='bottom',horizontalalignment='right')
#plt.text(4.e-2,5.e-13,'$p=6$',verticalalignment='top',horizontalalignment='left')
 
plt.tight_layout()
plt.savefig('./figures/prob1_convergence.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
plt.rc("font", size=16)
 
names = ['DeC%dGL'%(order) for order in range(2,10)]
dts = np.logspace(-2.5,0.,6)
u0 = np.array([1.,0.])
t_final = 10.
 
fig=plt.figure(figsize=(6,6))
ax = plt.subplot(111)
lines = []
 
for j, rkm in enumerate(decGLBRKs[:5]):
for relax in (False, True):
sols = []
for dt in dts:
 
tt, uu = RRK(rkm, dt,f_prob1, w0=u0, t_final=t_final, relaxation=relax, rescale_step=True, gammatol=1.0)
sols.append(uu[0,-1])
 
errors = []
ref = np.cos(t_final)
for i in range(len(sols)):
errors.append(abs(ref-sols[i]))
 
if relax:
plt.loglog(dts, errors,linestyle='dashed',color=colors[j],marker=markers[j],label=None);
else:
lines.append(plt.loglog(dts, errors,linestyle='-',color=colors[j],
marker=markers[j],label=names[j])[0])
 
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
 
# Put a legend to the right of the current axis
ax.legend(lines,names,loc='center left', bbox_to_anchor=(1, 0.5))
 
plt.xlabel('$\Delta t$')
plt.ylabel('Error');
plt.ylim([1.e-14, 10]);
 
plt.plot([1e-2,4e-2],[3.e-4, 4**2*3.e-4],'-ok',alpha=0.5)
plt.plot([1e-2,4e-2],[5.e-6, 4**3*5.e-6],'-ok',alpha=0.5)
plt.plot([1e-2,4e-2],[3.e-10, 4**4*3.e-10],'-ok',alpha=0.5)
plt.plot([1e-2,4.e-2],[5.e-13, 4**5*5.e-13],'-ok',alpha=0.5)
plt.plot([4e-2,16e-2],[5.e-13, 4**6*5.e-13],'-ok',alpha=0.5)
 
plt.text(1.e-2,3.e-4,'$p=2$',verticalalignment='bottom',horizontalalignment='right')
plt.text(1.e-2,5.e-6,'$p=3$',verticalalignment='top',horizontalalignment='left')
plt.text(3.e-2,6.e-9,'$p=4$',verticalalignment='bottom',horizontalalignment='left')
plt.text(1.e-2,5.e-13,'$p=5$',verticalalignment='bottom',horizontalalignment='right')
plt.text(4.e-2,5.e-13,'$p=6$',verticalalignment='top',horizontalalignment='left')
 
plt.tight_layout()
plt.savefig('./figures/prob1_DeCGLconvergence.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
dt = 0.5
order = 7
tt, uu = RRK(decEqRKs[order-2], dt, f_prob1, w0=u0, t_final=1000.,
relaxation=0, debug=False)
plt.plot(uu[0,:],uu[1,:],color='#56B4E9')
tt, uu = RRK(decEqRKs[order-2], dt, f_prob1, w0=u0, t_final=1000.,
relaxation=1, debug=False)
plt.plot(uu[0,:],uu[1,:],'-k')
plt.plot(uu[0,0],uu[1,0],'ok')
plt.axis('equal');
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
def makeplot(dt):
plt.figure(figsize=(8,8))
#dt = 1.5
unp1g, gam = RRK(dec3, dt, f_prob1, w0=u0, t_final=1000.,
relaxation=1, debug=False,one_step=True,gammatol=1.0)
du = (unp1g-u0)/gam
unp1 = u0+du
plt.plot([u0[0],unp1g[0]],[u0[1],unp1g[1]],'o-',color='#56B4E9')
plt.plot([u0[0],unp1[0]],[u0[1],unp1[1]],color='r')
 
t = np.linspace(0,2*np.pi,1000)
plt.plot(np.cos(t),np.sin(t),'-k',lw=0.5)
plt.plot(np.cos(dt),np.sin(dt),'ok')
plt.plot(np.cos(gam*dt),np.sin(gam*dt),'ob')
unorm = np.sqrt(unp1[0]**2+unp1[1]**2)
plt.plot(unp1[0]/unorm,unp1[1]/unorm,'og')
plt.axis('equal');
 
makeplot(1.5)
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
from ipywidgets import widgets, interact
 
interact(makeplot,dt=widgets.FloatSlider(min=0.1,max=1.5));
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
# Stability regions
 
%% Cell type:code id: tags:
 
``` python
plt.rc("font", size=22)
 
def relaxed_stability_polynomial_gammaval(rkm,gam):
P, Q = rkm.stability_function()
 
coeffs = P.coef
for i in range(len(P)):
coeffs[i] = gam*coeffs[i]
pg = np.poly1d(coeffs)
return pg
 
def plot_perturbed_stability_regions(rkm,gammas):
fig=plt.figure(figsize=(12,12))
q = np.poly1d([1])
for g in gammas:
p = relaxed_stability_polynomial_gammaval(rkm,g)
stability_function.plot_stability_region(p,q,filled=False,color='b',
fignum=fig.number,alpha=1.6-g**1.4,
bounds=[-3,0.6,-3.3,3.3])#,
#plot_axes=False);
p = relaxed_stability_polynomial_gammaval(rkm,1.)
stability_function.plot_stability_region(p,q,filled=False,color='k',bounds=[-5,5,-5,5],fignum=fig.number);
```
 
%% Cell type:code id: tags:
 
``` python
gammas = [0.7,0.8,0.9,1.0,1.1,1.2,1.3]
 
plot_perturbed_stability_regions(dec3,gammas)
plt.savefig('./figures/DeC3_stability_perturbed.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
plot_perturbed_stability_regions(dec4equi,gammas)
plt.savefig('./figures/DeC4equiv_stability_perturbed.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
plot_perturbed_stability_regions(dec4galo,gammas)
plt.savefig('./figures/DeC4galo_stability_perturbed.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
plot_perturbed_stability_regions(ssp3,gammas)
plt.savefig('./figures/RK3_stability_perturbed.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
plot_perturbed_stability_regions(rk4,gammas)
plt.savefig('./figures/RK4_stability_perturbed.pdf', bbox_inches="tight")
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
# Spectral advection
 
%% Cell type:code id: tags:
 
``` python
def RRK_with_spectrum(rkm, dt, f, w0=[1.,0], t_final=1., relaxation=True,
rescale_step=True, debug=False, gammatol=0.1):
w = np.array(w0)
t = 0
ww = np.zeros([len(w0),int((t_final-t)/dt*1.5)])
ww[:,0] = w.copy()
tt = [t]
ii = 0
b = rkm.b
s = len(rkm)
y = np.zeros((s,len(w0)))
max_gammam1 = 0.
gams = []
uhats = [np.abs(np.fft.fft(w0))]
while t < t_final:
if t + dt >= t_final:
dt = t_final - t
 
for i in range(s):
y[i,:] = w.copy()
for j in range(i):
y[i,:] += rkm.A[i,j]*dt*f(y[j,:])
 
F = np.array([f(y[i,:]) for i in range(s)])
 
if relaxation:
numer = 2*sum(b[i]*rkm.A[i,j]*np.dot(F[i],F[j]) \
for i in range(s) for j in range(s))
denom = sum(b[i]*b[j]*np.dot(F[i],F[j]) for i in range(s) for j in range(s))
if denom != 0:
gam = numer/denom
else:
gam = 1.
else: # Use standard RK method
gam = 1.
 
gm1 = np.abs(1.-gam)
max_gammam1 = max(max_gammam1,gm1)
if gm1 > gammatol:
print(gam)
raise Exception("The time step is probably too large.")
 
w = w + gam*dt*sum([b[j]*F[j] for j in range(s)])
if (t+dt < t_final) and rescale_step:
t += gam*dt
else:
t += dt
ii += 1
tt.append(t)
ww[:,ii] = w.copy()
gams.append(gam)
uhats.append(np.abs(np.fft.fft(w)))
if debug:
print(max_gammam1)
return tt, ww[:, :ii+1], np.array(gams), uhats
return tt, ww[:,:ii+1], uhats
```
 
%% Cell type:code id: tags:
 
``` python
plt.rc("font", size=20)
 
m=128
L = 2*np.pi
x = np.arange(-m/2,m/2)*(L/m)
xi = np.fft.fftfreq(m)*m*2*np.pi/L
dx = x[1]-x[0]
 
def f_advection(u):
uhat = np.fft.fft(u)
return -np.real(np.fft.ifft(1j*xi*uhat))
```
 
%% Cell type:code id: tags:
 
``` python
def plot_rel_amplification(rkname, IC='sech2', relaxation=True, mus=(0.2, 0.5, 0.9, 0.99), filename=None):
if IC == 'sech2':
u0 = 1./np.cosh(7.5*(x+1))**2
elif IC == 'whitenoise':
theta = np.random.rand(*xi.shape)
u0hat = np.exp(1j*theta)
u0 = np.real(np.fft.ifft(u0hat))
 
rkm = rk.loadRKM(rkname)
stabint = rkm.imaginary_stability_interval()
if stabint == 0: stabint = 1.
plt.figure(figsize=(6,5))
for mu in mus:
dt = mu * stabint * L/m/np.pi
tt, uu, uhats = RRK_with_spectrum(rkm.__num__(), dt, f_advection, w0=u0, t_final=1.,
relaxation=relaxation, debug=False, gammatol=1.0)
uhatdiff = np.ma.array(uhats[-1]-uhats[0])
uhatreldiff = uhatdiff/uhats[0]
plt.plot(np.sort(xi)[m//2:],uhatreldiff[np.argsort(xi)][m//2:],'-',lw=3,label='$\mu = ' + str(mu) + '$')
plt.xlabel('Wavenumber')
plt.ylabel('Relative amplification')
#plt.legend(['$\mu = ' + str(mu) + '$' for mu in np.array(mus)],loc=3);
plt.xlim(0,m/2-1)
plt.ylim(-1.2,0.8)
plt.tight_layout()
if filename:
plt.savefig(filename)
 
ax = plt.gca()
plt.figure()
handles, labels = ax.get_legend_handles_labels()
plt.figlegend(handles, labels, loc="center", ncol=5)
plt.savefig("./figures/waveeq_amp_legend.pdf", bbox_inches="tight");
```
 
%% Cell type:code id: tags:
 
``` python
plot_rel_amplification('RK44','whitenoise',relaxation=False,
filename='./figures/waveeq_RK44_standard.pdf')
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
plot_rel_amplification('RK44','whitenoise',relaxation=True,
filename='./figures/waveeq_RK44_whitenoise_relaxation.pdf')
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
plot_rel_amplification('RK44','sech2',relaxation=True,mus=(0.2, 0.5, 0.9, 0.99, 1.02),
filename='./figures/waveeq_RK44_sech2_relaxation.pdf')
```
 
%% Output
 
 
 
%% Cell type:code id: tags:
 
``` python
plot_rel_amplification('RK44','whitenoise',relaxation=True,
filename='./figures/waveeq_rk44_whitenoise_relaxation.pdf')
```
 
%% Output