Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Square Root Map for Falkner-Skan
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()
@jstults

This comment has been minimized.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.