#!/usr/bin/env python
# coding: utf-8

# # Advection-Diffusion Problems
# 

# * Consider the general PDE:
# 
# $$\phi_t + u\phi_x = \Gamma\phi_{xx},$$
# 
# Or
# 
# $$\frac{\partial\phi}{\partial t} + u\frac{\partial\phi}{\partial x} = \Gamma \frac{\partial^2\phi}{\partial x^2},$$
# 
# where $u$ is a velocity, $\phi$ is our scalar variable, and $\Gamma$ is a diffusion coefficient (like $\alpha$).
# 
# * The second term in the equation is a advection term. This is a hyperbolic-like term, a wave-like term.
# * The equation as a whole is parabolic.
# 
# 

# * Example. Let $\phi=\rho Y_i$, $\Gamma=D_i$ with constant $\rho$, $D_i$:
# $$\frac{\partial\rho Y_i}{\partial t} + \vec{v}\cdot\nabla(\rho Y_i) = \rho D_i\nabla^2Y_i.$$
#     * Unsteady species transport equation.

# ## FTCS method
# 
# $$\frac{\partial\phi}{\partial t} + u\frac{\partial\phi}{\partial x} = \Gamma \frac{\partial^2\phi}{\partial x^2},$$
# 
# $$\phi_i^{n+1} = \phi_i^n - 
# \frac{1}{2}\underbrace{\frac{u\Delta t}{\Delta x}}_{c}(\phi_{i+1}^n - \phi_{i-1}^n) + 
# \underbrace{\frac{\Gamma\Delta t}{\Delta x^2}}_{d}(\phi_{i-1}^n - 2\phi_i^n + \phi_{i+1}^n).$$
# 
# * FTCS properties: 
#     * consistent, 
#     * conditionally stable, 
#     * $\mathcal{O}(\Delta t)$, 
#     * $\mathcal{O}(\Delta x^2)$.
#     
# 

# * A stability analysis yeilds the following constraints:
# $$d\le \frac{1}{2},$$
# $$c^2 \le 2d,$$
# Or,
# $$c^2\le2d\le 1.$$
# 
# * For $d=1/2$, $c<1$.
# * Both $c$ and $d$ are time ratios.
#     * $d$ is the ratio of our timestep size to the grid-scale diffusion time.
#         * $\tau_D = \Delta x^2/\Gamma.$
#         * This is the nominal time it takes to diffuse across a cell.
#     * $c$ is the ratio of our timestep size to the grid-scale advection time.
#         * $\tau_c = \Delta x/u.$
#         * This is the nominal time it takes to convect across a cell.
#     * **Don't take a longer timestep than the physical time of advection or diffusion.**
# * Note that $c$ and $d$ are related.
# 

# ## BTCS method
# * BTCS properties: 
#     * consistent, 
#     * unconditionally stable, 
#     * $\mathcal{O}(\Delta t)$, 
#     * $\mathcal{O}(\Delta x^2)$.

# ## Trouble with central differences
# 
# * Consider the steady state equation
# 
# $$u\frac{\partial\phi}{\partial x} = \Gamma \frac{\partial^2\phi}{\partial x^2}.$$
# 
# * Apply a finite volume solution by integrating over a control volume:
# 
# 

# 
# $$(u\phi)_e - (u\phi)_w = \Gamma\left(\frac{d\phi}{dx}\right)_e - \Gamma\left(\frac{d\phi}{dx}\right)_w.$$
# 
# $$\frac{1}{2}u_e(\phi_E + \phi_P) - \frac{1}{2}u_w(\phi_P + \phi_W) = \frac{\Gamma}{\Delta x}(\phi_E - \phi_P) - \frac{\Gamma}{\Delta x}(\phi_P - \phi_W).$$
# 

# 
# * For constant $\rho$, $u$ is constant. Also assume $\Gamma$ is constant. Rearrange the expression:
# 
# $$\phi_W\left(-\frac{1}{2}u - \frac{\Gamma}{\Delta x}\right) + \phi_P\left(\frac{\Gamma}{\Delta x}\right) + \phi_E\left(\frac{1}{2}u - \frac{\Gamma}{\Delta x}\right) = 0.$$
# 
# Divide through by $\frac{\Gamma}{\Delta x}$:
# 
# $$\phi_W\left(-\frac{Pe}{2} -1 \right) + \phi_P(2) + \phi_E\left(\frac{Pe}{2} - 1\right) = 0.$$
# 
# * Here, $Pe=u\Delta x/\Gamma$ is a grid Peclet number. If $\Gamma=\nu$, the kinematic viscosity, this would be a grid Reynolds number. 
#     * $Pe = c/d = \tau_D/\tau_c$. 
#     * $Pe$ is a advection rate over a diffusion rate.
# 

# 
# Solve the above equation for $\phi_P$ in terms of its neighbors:
# 
# $$\phi_P = \frac{1}{2}\left[\phi_W\left(\frac{Pe}{2} + 1\right) - \phi_E\left(\frac{Pe}{2}-1\right)\right].$$
# 
# * Suppose we have a decreasing profile where $\phi_W>\phi_E$.
#     * For a steady state problem without a source term, we expect $\phi_P$ to be between its neighbors $\phi_E$ and $\phi_W$. 
#     * But if our velocity and grid spacing are set so that $Pe > 2$, then the sign of the $\phi_E$ term in the above equations changes and $\phi_P > \phi_W$, contrary to our expectation.
#     * For $Pe>2$, we get unphysical behavior.
#     * For given $u$ and $\Gamma$, $Pe$ can be reduced by decreasing $\Delta x$.
# 

# In[1]:


import ipywidgets as wgt
import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')


def f(Pe):
    ϕw = 2
    ϕe = 1
    ϕp = 1/2*(ϕw*(Pe/2 + 1) - ϕe*(Pe/2 - 1))
    ϕ = np.array([ϕw,ϕp,ϕe])
    x = np.array([0, 1, 2])
    plt.rc("font", size=14)
    plt.plot(x,ϕ,'o-', markersize=10, lw=3)
    plt.ylim([0,3])
    plt.xlabel('x')
    plt.ylabel('ϕ')
    plt.grid()
    plt.show() 

wgt.interact(f, Pe=(-5,5,0.1));     # try adding string and list arguments...


# Following Patankar, we can write our above equation as
# 
# $$C_p \phi_P = C_E\phi_E + C_W\phi_W,$$
# 
# where $C$'s are coefficients of the $\phi$'s.
# * Here, $C_P=1$, $C_E = -Pe/4 + 1/2$, and $C_W = Pe/4 + 1/2.$
# * All coefficients in this form should have the same sign so that increases in the neighbors result in increases in $\phi_P$, and vice-versa.
# * If $Pe>2$, then $C_E$ is negitive while the other $C$'s are positive. 

# ## Remedy: upwinding
# 
# The problem above is that we applied a central difference approximation when evaluating the convective term. 
# * (In the FV approach, we applied a linear interpolation when evaluating the face properties, which is analogous to using a central difference.)
# 
# These approaches effectively give equal weighting to the neighbors when computing a property. 
# * **But advection is directional.**
#     * You can smell much better when something is upwind of you then when it is downwind of you. 
# * In evaluating face properties, don't weight each neighbor equally. 
# * Insteady, advection *blows the upstream value to the face.*
# 
# $$\phi_e = \phi_P,\phantom{xx}\mbox{if}\phantom{xx} u_e>0,$$
# $$\phi_e = \phi_E,\phantom{xx}\mbox{if}\phantom{xx} u_e<0,$$
# 
# $$\phi_w = \phi_W,\phantom{xx}\mbox{if}\phantom{xx} u_w>0,$$
# $$\phi_w = \phi_P,\phantom{xx}\mbox{if}\phantom{xx} u_w<0.$$
# 
# 

# There is no change to the diffusion term. Diffusion has no directionality. Use central differences.
# 
# * Now, take $u>0$ and solve our PDE:
# 
# $$u\phi_e - u\phi_w = \Gamma\left(\frac{d\phi}{dx}\right)_e - \Gamma\left(\frac{d\phi}{dx}\right)_w.$$
# 
# $$u\phi_P - u\phi_W = \frac{\Gamma}{\Delta x}(\phi_E - \phi_P) - \frac{\Gamma}{\Delta x}(\phi_P - \phi_W).$$
# 
# 

# 
# Solve for $\phi_P$:
# 
# $$\phi_P = 
# \phi_E\left[\frac{\frac{\Gamma}{\Delta x}}{u + 2\frac{\Gamma}{\Delta x}}\right] + 
# \phi_W\left[\frac{u+\frac{\Gamma}{\Delta x}}{u + 2\frac{\Gamma}{\Delta x}}\right],
# $$
# 
# Or, 
# 
# $$\phi_P = 
# \phi_E\left[\frac{1}{Pe+2}\right] + 
# \phi_W\left[\frac{Pe+1}{Pe+2}\right].
# $$
# 
# * Here, the coefficients are always positive for any value of $Pe$.
# 
# * Note, this method is $\mathcal{O}(\Delta t)$, but only $\mathcal{O}(\Delta x)$. The upwinding is first order spatially.
#     

# In[ ]: