fdadvect.py

# finite-difference implementation of FTCS for linear advection.
#
# We are solving a_t + u a_x = 0
#
# The FTCS discretization is: anew = aold + (C/2) (aold_{i+1} - aold_{i-1})
#
# where C is the CFL number
#
# M. Zingale (2013-03-12)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['mathtext.fontset'] = 'cm'
mpl.rcParams['mathtext.rm'] = 'serif'

class FDGrid(object):

    def __init__(self, nx, ng, xmin=0.0, xmax=1.0):

        self.xmin = xmin
        self.xmax = xmax
        self.ng = ng
        self.nx = nx

        # python is zero-based.  Make easy intergers to know where the
        # real data lives
        self.ilo = ng
        self.ihi = ng+nx-1

        # physical coords
        self.dx = (xmax - xmin)/(nx-1)
        self.x = xmin + (np.arange(nx+2*ng)-ng)*self.dx

        # storage for the solution
        self.a = np.zeros((nx+2*ng), dtype=np.float64)

    def scratch_array(self):
        """ return a scratch array dimensioned for our grid """
        return np.zeros((self.nx+2*self.ng), dtype=np.float64)

    def fill_BCs(self):
        """ fill the a single ghostcell with periodic boundary conditions """
        self.a[self.ilo-1] = self.a[self.ihi-1]
        self.a[self.ihi+1] = self.a[self.ilo+1]


def solve_advection(g, u, C, method="upwind", tmax_factor=1.0):

    # time info
    dt = C*g.dx/u
    t = 0.0
    tmax = tmax_factor*(g.xmax - g.xmin)/u


    # initialize the data -- tophat
    g.a[:] = 0.0
    g.a[np.logical_and(g.x >= 1./3., g.x <= 2./3.)] = 1.0

    g.ainit = g.a.copy()

    # evolution loop
    anew = g.scratch_array()

    while t < tmax:

        # fill the boundary conditions
        g.fill_BCs()

        # loop over zones: note since we are periodic and both endpoints
        # are on the computational domain boundary, we don't have to
        # update both g.ilo and g.ihi -- we could set them equal instead.
        # But this is more general
        for i in range(g.ilo, g.ihi+1):

            if method == "upwind":
                anew[i] = g.a[i] - C*(g.a[i] - g.a[i-1])
            elif method == "FTCS":
                anew[i] = g.a[i] - 0.5*C*(g.a[i+1] - g.a[i-1])
            else:
                sys.exit("invalid method")

        # store the updated solution
        g.a[:] = anew[:]
        t += dt


# create the grid
nx = 65
ng = 1
g = FDGrid(nx, ng)


# define the CFL and speed
Clist = [0.1, 0.5, 0.9]
u = 1.0

# upwind
for n, C in enumerate(Clist):
    solve_advection(g, u, C)

    if n == 0:
        plt.plot(g.x[g.ilo:g.ihi+1], g.ainit[g.ilo:g.ihi+1], ls=":", label="exact")

    plt.plot(g.x[g.ilo:g.ihi+1], g.a[g.ilo:g.ihi+1], 
             label=r"$C = {}$".format(C))

plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$a$", fontsize=16)

plt.legend(frameon=False, loc="best")

plt.tight_layout()

plt.savefig("fdadvect-upwind.pdf")



plt.clf()

# upwind
for C in Clist:
    plt.clf()
    solve_advection(g, u, C, method="FTCS", tmax_factor=0.1)

    plt.plot(g.x[g.ilo:g.ihi+1], g.ainit[g.ilo:g.ihi+1], ls=":", label="exact")

    plt.plot(g.x[g.ilo:g.ihi+1], g.a[g.ilo:g.ihi+1], 
             label=r"$C = {}$".format(C))

    plt.xlabel(r"$x$", fontsize=16)
    plt.ylabel(r"$a$", fontsize=16)

    plt.legend(frameon=False, loc="best")

    plt.tight_layout()

    plt.savefig("fdadvect-FTCS-C{}.pdf".format(C))