Hyperbolic tangent point distribution for Falkner-Skan

arctanh_fs.py

import scipy as sp
from scipy.integrate import ode 
from scipy.optimize import fsolve 
from matplotlib import rc
rc('text', usetex=True)
tick_size = 'large'
label_size = 'xx-large'
rc('xtick', labelsize=8)
rc('ytick', labelsize=8)
rc('legend', fontsize=12)
rc('axes', labelsize=12)
fig_width_pt = 469.75502 
inches_per_pt = 1.0/72.27
golden_mean = (sp.sqrt(5.0)-1.0)/2.0
fig_width = fig_width_pt *inches_per_pt 
fig_height = fig_width * golden_mean
fig_size = [fig_width, fig_height] 
rc('figure', figsize=fig_size)

from matplotlib import pylab as plt

from fs import *

def int_tanh_fs(beta_0, beta, alpha, eta_inf):
    args = sp.array([beta_0, beta, eta_inf])    
    y0 = sp.array([0.0, 0.0, alpha])
    r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-14, atol=1e-14, method='adams', with_jacobian=True, order=12)
    t0 = 0.0 
    t1 = 1.0 
    # instead of using a uniform point distribution, use a tanh
    # mapping, which clusters the points near the wall where the
    # variation is largest
    nt = 20
    t = sp.arctanh(sp.linspace(t0,t1,nt,endpoint=False)) 
    dt = sp.diff(t) 
    r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
    for i in xrange(nt):
        r.integrate(r.t+dt[i])
    return(r)

def objective_func(x, beta_0, beta):
    alpha = x[0] 
    eta_inf = x[1]
    r = int_fs(beta_0, beta, alpha, eta_inf)
    return(sp.array([(r.y[1]-1.0).real, (r.y[2]).real]))

def extrap(y):
    return(2*y[1]-y[0])
def extrap2(y):
    return(-(4*y[1]-y[0]-5*y[2])/2.0)

foo = sp.genfromtxt("foo.txt")
beta = sp.genfromtxt("beta.txt")
#alpha = sp.genfromtxt("alpha.txt")
#eta_inf = sp.genfromtxt("eta_inf.txt")
beta_0 = 1.0

r = ode(fs.rates, fs.jac).set_integrator('zvode', rtol=1e-14, atol=1e-14, method='adams', with_jacobian=True, order=12)
nt = sp.logspace(4,8,5,base=2) 
u1 = []
t = []
dt = []
for j,nnt in enumerate(nt):
    u1.append(sp.zeros(nnt))
    u1[j][0] = 0.0 
    y0 = sp.array([0.0,0.0,foo[0,0]])
    t0 = 0.0 
    t1 = 1.0 
    # instead of using a uniform point distribution, use a tanh
    # mapping, which clusters the points near the wall where the
    # variation is largest
    t.append(sp.arctanh(sp.linspace(t0,t1,nnt,endpoint=False)))
    dt.append(sp.diff(t[j])) 
    args = sp.array([beta_0, beta[0], foo[0,1]])    
    r.set_initial_value(y0, 0.0).set_f_params(args).set_jac_params(args)
    for i in xrange(dt[j].shape[0]):
        r.integrate(r.t+dt[j][i])
        u1[j][i+1] = r.y[1].real  

plt.figure()
for i,u in enumerate(u1):
    plt.plot(u, foo[0,1]*t[i], label=r'$nt=%g$'%nt[i])
plt.legend(loc=0)
plt.xlabel(r"$f'$")
plt.ylabel(r'$\eta$')

plt.figure() # wall close-up
for i,u in enumerate(u1):
    plt.plot(u, foo[0,1]*t[i], '-x', label=r'$nt=%g$'%nt[i])
plt.legend(loc=0)
plt.xlabel(r"$f'$")
plt.ylabel(r'$\eta$')
xmin,xmax,ymin,ymax = plt.axis()
plt.axis([xmin,0.017,ymin,0.45])
plt.show()

Not really a grid convergence study since the error control on the integrator is cranked down so tight.
Close-up of near-wall region