Total Variation Diminishing (TVD) schemes are flux limiters for the convective term that keep an advected profile bounded and non-oscillatory near a discontinuity, without the excessive numerical diffusion of first-order upwind. This notebook compares several limiters (upwind, Min-mod, Osher, CLAM, MUSCL, van Leer) on the same 1D linear advection test.
Limiter functions (NVD diagrams)¶
Each limiter is a function of the normalized upstream value that returns the normalized downwind-face value (plotted here as ), for a fixed value of the local Courant-like ratio. This is the Normalized Variable Diagram (NVD) representation of a limiter, independent of any PDE; it shows how aggressively each scheme interpolates towards the downwind cell.
import numpy as np
import matplotlib.pyplot as plt
from pymrm import upwind, minmod, osher, clam, muscl, vanleer
limiters= [
(upwind, "upwind"),
(minmod, "Min-mod"),
(osher, "Osher"),
(clam, "CLAM"),
(muscl, "MUSCL"),
(vanleer, "van Leer")
]
c=np.linspace(-1,2,99)
for func, name in limiters:
plt.plot(c, func(c, 0.5, 0.75), label=name)
plt.title('NVD diagrams limiter functions')
plt.xlabel(r'$\tilde{c}_C$')
plt.ylabel(r'$\tilde{c}_d - \tilde{c}_C$')
plt.legend()
plt.show()
Explicit time stepping: step advection¶
Pure advection at velocity of a step profile: a Danckwerts-type inlet at fixes , and a zero-gradient outlet at lets the profile leave without reflection. Starting from , each limiter is advanced with explicit (Euler-forward) time stepping and the final profile is plotted.
import math
import numpy as np
import matplotlib.pyplot as plt
from pymrm import interp_cntr_to_stagg_tvd
from pymrm import upwind, minmod, osher, clam, muscl, vanleer
limiters= [
(upwind, "upwind"),
(minmod, "Min-mod"),
(osher, "Osher"),
(clam, "CLAM"),
(muscl, "MUSCL"),
(vanleer, "van Leer")
]
L = 10;
num_x = 20;
v = 1
Co = 0.25;
dt = Co*(L/num_x)/v
num_time_steps = math.floor(0.75*num_x/Co)
x_f = np.linspace(0,L,num_x+1)
x_c = 0.5*(x_f[1:]+x_f[:-1])
bc = ({'a':0, 'b':1, 'd':1},{'a':1, 'b':0, 'd':0})
for func, name in limiters:
c = np.zeros(num_x)
for i in range(num_time_steps):
c_f, dc_f = interp_cntr_to_stagg_tvd(c, x_f, x_c, bc, v, tvd_limiter=func, axis=0)
c += Co*(c_f[:-1]-c_f[1:])
plt.plot(x_c, c, label=name)
plt.title(fr"Euler-forward explicit convection, $Co=$ {Co}")
plt.xlabel(r'position, $x$')
plt.ylabel(r'concentration, $c$')
plt.legend()
plt.show()
Implicit time stepping with deferred correction¶
The same test, now advanced with implicit (Euler-backward) time stepping. The Jacobian is built once from the first-order upwind convection operator; the TVD limiter only enters as an explicit correction term recomputed each inner iteration (deferred correction), so the constant operators and their LU factorization are reused for every limiter and time step.
import math
import numpy as np
from scipy.sparse import eye_array
from scipy.sparse import linalg as sla
from pymrm import construct_convflux_upwind, construct_div, interp_cntr_to_stagg_tvd
from pymrm import upwind, minmod, osher, clam, muscl, vanleer
limiters= [
(upwind, "upwind"),
(minmod, "Min-mod"),
(osher, "Osher"),
(clam, "CLAM"),
(muscl, "MUSCL"),
(vanleer, "van Leer")
]
L = 10
num_x = 20
v = 1
bc = ({'a':0, 'b':1, 'd':1},{'a':1, 'b':0, 'd':0})
Co = 0.25
dt = Co*(L/num_x)/v
num_time_steps = math.floor(0.75*num_x/Co)
x_f = np.linspace(0,L,num_x+1)
x_c = 0.5*(x_f[1:]+x_f[:-1])
# jac_mat = (1/dt*I + div_mat @ conv_mat)
# g = (c-c_old)/dt + div_mat @ (conv_mat @ c + conv_bc + v*dc_f)
# = jac_mat @ c - c_old/dt + div_mat @ (conv_bc + v*dc_f)
num_inner_iter = 1
conv_mat, conv_bc = construct_convflux_upwind(num_x, x_f, x_c, bc, v, axis=0)
div_mat = construct_div(num_x, x_f, nu=0, axis=0)
jac_mat = (1/dt)*eye_array(num_x, format='csc') + div_mat @ conv_mat
jac_lu = sla.splu(jac_mat)
for func, name in limiters:
c = np.zeros(num_x)
for i in range(num_time_steps):
c_old = c.copy()
for j in range(num_inner_iter):
c_f, dc_f = interp_cntr_to_stagg_tvd(c, x_f, x_c, bc, v, tvd_limiter=func, axis=0)
g = jac_mat @ c.reshape(-1,1) - c_old.reshape(-1,1)/dt + div_mat @ (conv_bc + v*dc_f.reshape(-1,1))
c -= jac_lu.solve(g).reshape(c.shape)
plt.plot(x_c, c, label=name)
plt.title(fr"Euler-backward implicit convection with deferred correction, $Co=$ {Co}")
plt.xlabel(r'position, $x$')
plt.ylabel(r'concentration, $c$')
plt.legend()
plt.show()
Checking the result¶
TVD limiters are designed to be boundedness-preserving (Total Variation Diminishing): with an inlet value of 1, an initial value of 0, and no source term, the concentration must stay within at every cell and every time step, for every limiter. A violation of this bound would mean the limiter, or its implementation, is not actually TVD.
import math
import numpy as np
from pymrm import interp_cntr_to_stagg_tvd
from pymrm import upwind, minmod, osher, clam, muscl, vanleer
limiters= [
(upwind, "upwind"),
(minmod, "Min-mod"),
(osher, "Osher"),
(clam, "CLAM"),
(muscl, "MUSCL"),
(vanleer, "van Leer")
]
L = 10;
num_x = 20;
v = 1
Co = 0.25;
dt = Co*(L/num_x)/v
num_time_steps = math.floor(0.75*num_x/Co)
x_f = np.linspace(0,L,num_x+1)
x_c = 0.5*(x_f[1:]+x_f[:-1])
bc = ({'a':0, 'b':1, 'd':1},{'a':1, 'b':0, 'd':0})
tol = 1e-10
print(f"{'limiter':<10s} {'min(c)':>12s} {'max(c)':>12s} interpretation")
for func, name in limiters:
c_chk = np.zeros(num_x)
c_min, c_max = c_chk.min(), c_chk.max()
for i in range(num_time_steps):
c_f, dc_f = interp_cntr_to_stagg_tvd(c_chk, x_f, x_c, bc, v, tvd_limiter=func, axis=0)
c_chk += Co*(c_f[:-1]-c_f[1:])
c_min = min(c_min, c_chk.min())
c_max = max(c_max, c_chk.max())
bounded = c_min >= -tol and c_max <= 1 + tol
verdict = "bounded" if bounded else "NOT bounded"
print(f"{name:<10s} {c_min:12.4e} {c_max:12.4e} {verdict}")
limiter min(c) max(c) interpretation
upwind 0.0000e+00 1.0000e+00 bounded
Min-mod 0.0000e+00 1.0000e+00 bounded
Osher -2.4449e-17 1.0000e+00 bounded
CLAM -2.4449e-17 1.0000e+00 bounded
MUSCL -2.3235e-17 1.0000e+00 bounded
van Leer -8.0861e-26 1.0000e+00 bounded