Square Root Map for Falkner-Skan

## integrate_fs3.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=14)
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_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
    nt = 200
    dt = t1 / float(nt)
    r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
    for i in xrange(nt):
        r.integrate(r.t+dt)
    return(r)

def int_tanh_fs(beta_0, beta, alpha, eta_inf, nt):
    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 = 64
    t = sp.arctanh(sp.linspace(t0,t1,nt,endpoint=False))
#    t = t/t.max()
    dt = sp.diff(t)
    r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
    for i in xrange(nt-1):
        r.integrate(r.t+dt[i])
    return(r)

def int_sqrt_fs(beta_0, beta, alpha, eta_inf, nt):
    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=2e-16, atol=2e-16, method='adams', with_jacobian=True, order=12)
    t0 = 0.0
    t1 = 1.0
    # instead of using a uniform point distribution, use a square root
    # mapping, which clusters the points near the wall where the
    # variation is largest
#    nt = 64
    x = sp.linspace(t0,t1,nt,endpoint=False)
    t = y_of_x(eta_inf,x)
#    t = t/t.max()
    dt = sp.diff(t)
    r.set_initial_value(y0, t0).set_f_params(args).set_jac_params(args)
    for i in xrange(int(nt)-1):
        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 objective_func2(x, beta_0, beta, nt):
    alpha = x[0]
    eta_inf = 1.0 # x[1]
    r = int_tanh_fs(beta_0, beta, alpha, eta_inf, nt)
    return(sp.array([(r.y[1]-1.0).real, (r.y[2]).real]))

def objective_func3(x, beta_0, beta, nt):
    alpha = x[0]
    eta_inf = 1.0 # x[1]
    r = int_sqrt_fs(beta_0, beta, alpha, eta_inf, nt)
    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)
# sqrt mapping
def y_of_x(L,x):
    return(L*x / sp.sqrt(1 - x**2))

eta_inf = 1.0
alpha = 1e-2
beta_0 = 1.0
x = sp.array([alpha, eta_inf])
nnt = sp.unique(sp.logspace(3,12,73,base=2).astype(int))
foo = sp.zeros((nnt.shape[0],2))
deta0 = sp.zeros(nnt.shape[0])
for i,nt in enumerate(nnt):
    args = (beta_0, 0.0, int(nt))
    foo[i] = fsolve(objective_func3, x, args=args, xtol=2e-16)
    if(i==0):
        x[0] = alpha
    elif(i==1):
        x[0] = extrap(foo[i-1:i+1,0])
    else:
        x[0] = extrap2(foo[i-2:i+1,0])
    x[1] = eta_inf
    y = y_of_x(eta_inf, sp.linspace(0,1,nt,endpoint=False))
    foo[i,1] = y.max()
    deta0[i] = y[1] - y[0] # initial point spacing off the wall
err = foo[:,0] - foo[foo.shape[0]-1,0]

plt.figure()
plt.loglog(nnt[0:nnt.shape[0]-1], abs(err[0:nnt.shape[0]-1]),label=r'$\eta=\frac{x}{\sqrt{1 - x^2}}$, $x \in [0,1)$')
plt.xlabel(r'$n$')
plt.ylabel(r'$\| \alpha_n - \alpha \|$')
plt.legend(loc=0)
plt.savefig("sqrt_map_cnvrg.png")
plt.figure()
plt.plot(nnt, foo[:,1],label=r'$\eta=\frac{x}{\sqrt{1 - x^2}}$, $x \in [0,1)$')
plt.xlabel(r'$n$')
plt.ylabel(r'$\eta_{\infty}$')
plt.legend(loc=0)
plt.savefig("eta_inf.png")
plt.figure()
plt.semilogy(nnt, deta0,label=r'$\eta=\frac{x}{\sqrt{1 - x^2}}$, $x \in [0,1)$')
plt.xlabel(r'$n$')
plt.ylabel(r'$\Delta \eta_{wall}$')
plt.legend(loc=0)
plt.savefig("deta0.png")
plt.show()