Acoustics with Godunov's method

R.J. LeVeque, AMath 574, Winter Quarter 2023

A simple implementation of Godunov's method for constant coefficient acoustics.

Modified to include dt in function calls, and to better illustrate the solution after 1 or more steps.

%matplotlib inline

from pylab import *

from IPython.display import HTML

try:
    from clawpack.visclaw import animation_tools
except:
    print("Failed to load animation_tools from Clawpack")

Define the domain and material parameters for acoustics:

# domain:
xlower = 0.
xupper = 1.

# material parameters:
K0 = 2
rho0 = 2
c0 = sqrt(K0/rho0)
print('sound speed c0 = %.3f' % c0)
Z0 = sqrt(K0*rho0)
print('impedance Z0 = %.3f' % Z0)

sound speed c0 = 1.000
impedance Z0 = 2.000

Function to take a single time step:

For simplicity, this function has only two inputs, the arrays of Pn $=P^n$ and Un $=U^n$ containing the pressure and velocity at some time t_n.

It returns the updated arrays Pnp $= P^{n+1}$ and Unp $= U^{n+1}$ at the end of the time step.

This function uses values set above for the material parameters set above.
It should work regardless of the number of grid cells in the arrays (provided Pn and Un have the same length), and now take dt as a parameter. So you only need to re-execute if you change the material parameters.

def Godunov_step(Pn,Un,dt):
    # initialize new cell averages at n+1 to old ones at n:
    # important to make a copy so we can modify without change old values!
    Pnp = Pn.copy()
    Unp = Un.copy()
    
    # fill ghost cells for periodic BCs:
    # fill at t_n so these values can be used in computing updates
    Pn[0] = Pn[-2]
    Pn[-1] = Pn[1]
    Un[0] = Un[-2]
    Un[-1] = Un[1]
    
    # loop over interfaces where Riemann problem will be solved:
    for i in range(1,len(x)):
        
        # solve Riemann problem between i-1 and i:
        
        # left and right states (using old values at t_n):
        PL = Pn[i-1]; PR = Pn[i]
        UL = Un[i-1]; UR = Un[i]
        
        # coefficients of eigenvectors in W^p = alpha^p * r^p:
        # from (3.31) in FVMHP
        alpha1 = (-(PR-PL) + Z0*(UR-UL)) / (2*Z0)
        alpha2 = ((PR-PL) + Z0*(UR-UL)) / (2*Z0)
        
        # wave speeds (eigenvalues):
        s1 = -c0
        s2 = c0
        
        # left-going wave W1 updates cell i-1:

        Pnp[i-1] = Pnp[i-1] - dt/dx * s1 * alpha1 * (-Z0)
        Unp[i-1] = Unp[i-1] - dt/dx * s1 * alpha1 * 1.
        
        # right-going wave W2 updates cell i:
        Pnp[i] = Pnp[i] - dt/dx * s2 * alpha2 * Z0
        Unp[i] = Unp[i] - dt/dx * s2 * alpha2 * 1.
        
    return Pnp, Unp

Function to plot a numerical solution at one time:

def plotQ(Pn,Un,tn):
    subplot(2,1,1)
    plot(x[1:-1], Pn[1:-1], 'bo-')
    xlim(xlower,xupper)
    title('pressure at time %.3f' % tn)
    
    subplot(2,1,2)
    plot(x[1:-1], Un[1:-1], 'bo-')
    xlim(xlower,xupper)
    title('velocity at time %.3f' % tn)
    tight_layout() 

Time stepping:

Executing the next cell repeatedly will take a time step and plot the new solution

def time_stepper(t0, P0, U0, dt, nsteps, one_step=Godunov_step):
    """
    Take nsteps with time step dt, starting with initial data (P0,U0) at time t0.
    To take a single step use the method specified by one_step, which
    defaults to Godunov_step since that is the only method currently defined.
    But you might want to define a new method for comparison.
    """
    Pn = P0.copy()
    Un = U0.copy()
    for n in range(nsteps):
        Pn,Un = one_step(Pn, Un, dt)
    return Pn,Un

Perform some numerical experiments

Now that things are set up we can do some experiments with various choices of the grid resolution, Courant number, initial conditions, etc.

Define the numerical grid

We define a finite volume grid with num_cells interior points, indexed 1,2,...,num_cells. The cell centers start at dx/2. There is also one "ghost cell" on either side. This works great in Python since array indexing starts at 0, but in Matlab with 1-based indexing the ghost cell would have index 1. Also if more than 1 ghost cell is used, the indexing will change.

num_cells = 50
dx = (xupper - xlower)/num_cells

# cell centers, including one ghost cell on either side:
x = arange(xlower-dx/2, xupper+dx, dx)

print('including 2 ghost cells, the grid has %i cells' % len(x))

including 2 ghost cells, the grid has 52 cells

Initial conditions:

The arrays P0 and U0 should be numpy arrays of length num_cells + 2 with dtype = float. Make sure they aren't integer arrays!

#P0 = exp(-150*(x-0.5)**2)
P0 = where(logical_and(x>0.4,x<0.6), 1., 0.)

U0 = zeros(x.shape)
#U0 = -P0 / Z0  # only left-going wave
#U0 = P0 / Z0  # only right-going wave

t0 = 0.
fig = figure(figsize=(5,5))
plotQ(P0,U0,t0)

Take several time steps and plot final solution:

dt = 0.018
cfl = c0*dt/dx
nsteps = 20
tn = t0 + nsteps*dt

print('Using dt = %.4f, cfl = %.2f, taking %i steps to time %.3f' \
      % (dt,cfl,nsteps,tn))

Pn,Un = time_stepper(t0, P0, U0, dt, nsteps, Godunov_step)
fig = figure(figsize=(5,5))
plotQ(Pn,Un,tn)

Using dt = 0.0180, cfl = 0.90, taking 20 steps to time 0.360

Second experiment - refine the grid

Try the same test but with 100 cells instead of 50, with the time step also reduced by a factor of 2, and now taking 40 time steps to reach the same time as with 20 steps before.

Now all the setup and plotting has been combined into a single cell...

num_cells = 100
dx = (xupper - xlower)/num_cells

# cell centers, including one ghost cell on either side:
x = arange(xlower-dx/2, xupper+dx, dx)

print('including 2 ghost cells, the grid has %i cells' % len(x))

#P0 = exp(-150*(x-0.5)**2)
P0 = where(logical_and(x>0.4,x<0.6), 1., 0.)

U0 = zeros(x.shape)
#U0 = -P0 / Z0  # only left-going wave
#U0 = P0 / Z0  # only right-going wave

t0 = 0.
fig = figure(figsize=(5,5))
plotQ(P0,U0,t0)

dt = 0.009
cfl = c0*dt/dx
nsteps = 40
tn = t0 + nsteps*dt
print('Using dt = %.4f, cfl = %.2f, taking %i steps to time %.3f' \
      % (dt,cfl,nsteps,tn))

Pn,Un = time_stepper(t0, P0, U0, dt, nsteps, Godunov_step)
fig = figure(figsize=(5,5))
plotQ(Pn,Un,tn)

including 2 ghost cells, the grid has 102 cells
Using dt = 0.0090, cfl = 0.90, taking 40 steps to time 0.360

Make an animation:

This will only work if you were able to import animation_tools.

# Using grid points, dt, cfl, etc. from above.

# initial conditions from above:
tn = 0.
Pn = P0.copy()
Un = U0.copy()

nsteps = 40
nplot = 4  # how frequently to produce a plot
figs = []  # to accumulate figures for animation

for n in range(0,nsteps):
    tn = n*dt
    if mod(n,nplot)==0:
        fig = figure(figsize=(5,5))
        plotQ(Pn,Un,tn)
        figs.append(fig)
        close(fig)
    if n < nsteps-1:
        # take the next step
        Pn, Un = Godunov_step(Pn,Un,dt)


anim = animation_tools.animate_figs(figs, figsize=(5,5))
HTML(anim.to_jshtml())

Once Loop Reflect