fipy test of convection terms

test1d_convection.py

import fipy as fp
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

L = 1.0  # domain size
nx = 400  # discretization
dx = L/nx  # h

v = 0.1  # sampled coefficient

# mesh
m = fp.Grid1D(dx=dx, nx=nx)

# analytical solution at cell centers
axis = 0
x = np.array(m.cellCenters[axis],copy=True)
u = (np.exp(v*x) - np.exp(v))/(1.0 - np.exp(v))

fig2 = plt.figure(2,figsize=(8,5))
ax1 = fig2.add_subplot(121)
ax2 = fig2.add_subplot(122)

# available fipy convection terms
# http://www.ctcms.nist.gov/fipy/documentation/numerical/scheme.html
# ----------------------------------------
# ExponentialConvectionTerm <-- analytical solution for this problem
# CentralDifferenceConvectionTerm
# HybridConvectionTerm
# PowerLawConvectionTerm
# UpwindConvectionTerm


terms = [fp.ExponentialConvectionTerm, fp.CentralDifferenceConvectionTerm,
         fp.HybridConvectionTerm, fp.PowerLawConvectionTerm,
         fp.UpwindConvectionTerm]

terms = [fp.CentralDifferenceConvectionTerm]

colors = ['red','green',
          'blue','cyan',
          'magenta','black']

for color,ConvectionTerm in zip(colors,terms):

    print ConvectionTerm.__name__

    # variable (initial guess)
    U = fp.CellVariable(name='U', mesh=m)

    # boundary conditions
    U.constrain(1.0, m.facesLeft)  # u(0)=1
    U.constrain(0.0, m.facesRight) # u(1)=0
    # U.faceGrad.constrain(0.0, m.facesLeft)
    # U.faceGrad.constrain(0.1, m.facesRight)

    # governing equation
    eq = fp.DiffusionTerm(coeff=1.0,var=U) == ConvectionTerm(coeff=[v,],var=U)

    # steady-state solve
    eq.solve(var=U, solver=fp.LinearLUSolver())

    ax1.plot(x, U.value, linestyle='-', color=color)
    ax1.plot(x, U.value, linestyle='None', color=color, marker='.')

    ax2.semilogy(x, np.abs(u-U.value), linestyle='-', color=color,
                 label=ConvectionTerm.__name__.rstrip('ConvectionTerm'))
    ax2.semilogy(x, np.abs(u-U.value), linestyle='None', color=color, marker='.')

ax1.plot(x,u,'r-',label='analytical')
ax1.set_title('comparison')
ax1.set_xlabel('$x$')
ax1.set_ylabel('$u(v=%.1f,x)$' % v)

ax2.legend(loc=0,fontsize=8)
ax2.set_title('analytical - N(%i)' % nx)
ax2.set_xlabel('$x$')
ax2.set_ylabel('$|u(v,x)-U|$')

plt.tight_layout()
plt.savefig('comparison_N%i.png' % nx)
plt.close(2)