CFD Code Snippets—3: First-order wave equation, first-order upwind scheme

The first-order (i.e. one-way) wave equation in 1D is given by:

\dfrac{\partial T}{\partial t} = c \dfrac{\partial T}{\partial x}.

where T is the transverse displacement now, not temperature, and c is a constant.

The equation describes a pulse moving at a constant speed c to the right. If c < 0, the pulse moves to the left. Gilbert Strang [^] describes this equation with terms like (momentary) “flash of light” and “wall of water.”

It is with this problem that Prof. Lorena Barba chooses to begin [^] her 12 Steps program [^]. However, I am not so sure whether it is better to begin the course with convection rather than with diffusion.

The differential terms in convection are all of the first-order whereas those in the diffusion are mixed: first-order in time and second-order in space. It would thus appear that convection is a simpler phenomena than diffusion. In a certain sense, it is.

However, it so happens that for a course on CFD, the linear convection is only a stepping stone to the nonlinear convection. Seen from this viewpoint, it’s better to first cover diffusion, and only then go over to convection. This way, it becomes easier for the student to compare and contrast the nonlinear convection from the linear convection.

Currently, one way or the other, I am constrained to teach only according to the SPPU syllabus. So, one way or the other, the first-order wave equation could have come only after the diffusion equation.

No, don’t give too much credit to  the S. P. Pune University. The sequence wouldn’t be a conscious choice of the syllabus framers. The most probable scenario is that they must have just arranged and rearranged the material until they were satisfied that every committee-member’s suggestion was taken into consideration without any consideration of merit, and that all the six units carried more or less the same amount of material. Or, less likely, until the point that they themselves got tired of arguing with each other. In any case, none could possibly have thought of ordering the material according to a proper hierarchy.

Anyway, here we go with the Python code. (I hurriedly rewrote it afresh after the recent HDD crash, and so, comments and all could be refined a bit further. As to bugs: You be the tester.)


#################################################################
'''
First-Order Wave Equation, using First-Order Upwind Scheme
Code written by and copyright (c)Ajit R. Jadhav.

Based on the code by Prof. Lorena Barba (12 Steps to NS)
For theory, also see notes/resources on the 'net, such as:
Gilbert Strang: http://ocw.mit.edu/courses/mathematics/18-086-mathematical-methods-for-engineers-ii-spring-2006/readings/am52.pdf
Olver: http://www.math.umn.edu/~olver/num_/lnp.pdf
Gardner: https://math.la.asu.edu/~gardner/wave.pdf
and others.

Version: 29-March-2016.
'''
#################################################################

import numpy
from matplotlib import pyplot as plt
import random

#################################################################
#  Program Data. All units are in SI
#

# Physical time over which the simulation is to be run, in s
totalTime = 2.0

# Domain size, in meters
domLen = 1.0

# Signal velocity of the pulse, m/s
c = 1.0

# Number of cells
nc = 50
# Number of nodes
nx = nc + 1
# Compute the width of each cell
dx = domLen / nc

# Vertical (i.e. transverse) wave displacements in m, to use
# in constructing the initial and boundary conditions
T_Level_1 = 0.0
T_Level_2 = 0.1

# Create the numpy arrays, initialized to zero, for holding
# the primary variable values (here, transverse displacements)
# The main array, holds the values after time-marching through one step
T = numpy.zeros(nx)
# A temporary array, holds the values at the previous time-step
Tn = numpy.zeros(nx)

# Compute the duration of each time-step, and the number of time-steps required
# Stability criterion: CFL \equiv c \dfrac{\Delta t}{\Delta x}
# CFL <= 1.0 for stability, but the value of 1.0 here is special!
# Try other values of CFL < 1.0, a few values > 1.0, and a value of exactly 1.0
# Explain the smoothening at CFL < 1.0 # Explain the instability at CFL > 1.0
# Explain the ``oddity'' of the result at CFL = 1.0
CFL = 0.25

print("Assumed stability criterion value: %lf" % CFL)
# Given the above CFL value, back-calculate the time-step duration
dt = CFL * dx / c
print("Duration of each time-step in s, after applying the stability criterion: %lf" % dt)
# Compute the total no. of time-steps (an integer) to fit within the given time-interval
nt = int(totalTime / dt)
print("Number of time-steps required to step through: %d" % nt)

#################################################################
#  Setting up of the initial condition:
#
#  Uncomment the required IC and keep others commented

#########################
# IC1: A single pulse in the middle
# TL = T_Level_1
# TR = T_Level_1
# T[ nx/2 ] = T_Level_2

#########################
# IC2: A step function in the middle
# TL = T_Level_1
# TR = T_Level_1
# T[ nx/3.0 : 2.0*nx/3.0  ] = T_Level_2

#########################
# IC3: Left half at tempLevel_1; right half at tempLevel_2
# TL = T_Level_1
# TR = T_Level_2
# T[0: nx / 2] = TL
# T[nx / 2 + 1: nx] = TR

#########################
# IC4: An arbitrary function
random.seed(10)
TL = T_Level_1
TR = T_Level_1
for i in range(nc/5,nc/5+5):
    T[i] = T_Level_2 * random.random()

print("Initial condition is:")
print(T)
#
# Setting up of the initial condition is over
#################################################################

# Create the x-axis vector, for plotting purposes
xAxis = numpy.linspace(0, domLen, nx)

# Total time elapsed---physical time
tt = 0.0
plt.title("1D First-Order Wave Equation using\nFirst-Order Upwind Scheme. CFL: %5.2lf. Time: %5.2lf s" % (CFL, tt))
plt.xlabel("Distance, in m")
plt.ylabel("Transverse Displacement, in m")
plt.xlim(0.0, domLen)
plt.ylim(T_Level_1, T_Level_2 * 1.5)
plt.plot(xAxis, T, 'r')

#################################################################
#  Time-stepping, using upwind scheme
#

for n in range(nt):
    Tn = T.copy()
    # The trick on the following line is specific to Python's negative indexing
    for i in range(0, nx):
        # The following commented line would be the usual way of writing the loop
        # for i in range(1,nx):
        # The following step is where the first-order upwind scheme computation is done
        # Modify this one line to implement other schemes like Lax-Wendroff etc.
        T[i] = Tn[i] - CFL * (Tn[i] - Tn[i - 1])

    tt = dt * (n + 1)
    plt.clf()
    plt.title("1D First-Order Wave Equation using\nFirst-Order Upwind Scheme. CFL: %5.2lf. Time: %5.2lf s" % (CFL, tt))
    plt.xlabel("Distance, in m")
    plt.ylabel("Transverse Displacement, in m")
    plt.xlim(0.0, domLen)
    plt.ylim(T_Level_1, T_Level_2 * 1.5)
    plt.plot(xAxis, T, 'k')

    plt.pause(0.1)

plt.plot(xAxis, T);
plt.show()

print("done")


Some things to try:

  1. Straight-forward: As mentioned in the code, change the CFL to various values all of which are strictly less than 1.0. Then, change CFL to a few values greater than 1.0. Draw conclusions about stability.
  2. Not so straight-forward: Then, also set the CFL exactly to 1.0, and see what happens. Did you expect the outcome? Why or why not?
  3. Programming-related: Implement the Lax-Wendroff/MacCormack scheme.
  4. The smoothing that occurs at CFL values < 1.0 is due to a phenomenon called false diffusion. Explain the similarity and differences of the false diffusion from the true diffusion (say as modeled in the steps 1 and 2 of this series).

A Song I Like:

(Marathi) “dis chaar jhaale mann…”
Music: Ashok Patki
Lyrics: Soumitra
Singer: Sadhana Sargam


[Minor editing may happen…]

[E&OE]

Advertisements