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

# In[1]:


# Custom styling, borrowed from https://github.com/ellisonbg/talk-2013-scipy

from IPython import display


# Make plots centered. H/T to http://stackoverflow.com/a/27168595/1068170
STYLE = """\
<style>
div.output_area, .ui-wrapper {
  margin-left: auto !important;
  margin-right: auto !important;
}
</style>
"""


SLIDES_STYLE = """
<style>

h1.bigtitle {
    margin: 4cm 1cm 4cm 1cm;
    font-size: 300%;
}

h1.title { 
    font-size: 250%;
}

.rendered_html h1 {
    margin: 0.25em 0em 0.5em;
    color: #015C9C;
    text-align: center;
    line-height: 1.2; 
    page-break-before: always;
}
</style>
"""

STYLE = STYLE + '\n' + SLIDES_STYLE

display.display(display.HTML(STYLE))


# $$u_t + \left[f(u)\right]_x = 0$$

# <h1 class="bigtitle">Finite Volume</h1>

# $$w(x, t) = \frac{1}{\Delta x} \int_{x - \Delta x / 2}^{x + \Delta x / 2} u(s, t) \, ds$$

# $$\begin{align*}
# w_t &= \frac{1}{\Delta x} \int_{x - \Delta x / 2}^{x + \Delta x / 2} u_t \, ds \\
# &= \frac{1}{\Delta x} \int_{x - \Delta x / 2}^{x + \Delta x / 2} -\left[f(u)\right]_s \, ds \\
# &= -\frac{f\left(u\left(x + \Delta x / 2, t\right)\right) - f\left(u\left(x - \Delta x / 2, t\right)\right)}{\Delta x}
# \end{align*}$$

# $$\begin{align*}
# u_j^n &= u\left(j \Delta x, n \Delta t\right) \\
# w_j^n &= w\left(j \Delta x, n \Delta t\right)
# \end{align*}$$

# $$\frac{d}{dt} w_j^n = -\frac{1}{\Delta x}\left(f\left(u_{j + \frac{1}{2}}\right) - f\left(u_{j - \frac{1}{2}}\right)\right)$$

# $$\left\{w_j^n\right\}_j \overset{?}{\longrightarrow} \left\{u_{j + \frac{1}{2}}^n\right\}_j$$

# > WENO interpolation ... [is] used ... to transfer information from one domain to another in a high order, nonoscillatory fashion
# 
# ![WENO Transfer Quote](weno_transfer_quote.png "WENO Paper Quote")

# <h1 class="bigtitle">Finite Difference</h1>

# $$\frac{d}{dt} u_j^n = -\frac{1}{\Delta x}\left(\widehat{f}_{j + \frac{1}{2}} - \widehat{f}_{j - \frac{1}{2}}\right)$$

# Define $h(x)$ implicitly via
# $$\frac{1}{\Delta x} \int_{x - \Delta x / 2}^{x + \Delta x / 2} h(s) \, ds = f(u(x))$$

# $$\frac{1}{\Delta x} \int_{x - \Delta x / 2}^{x + \Delta x / 2} h(s) \, ds = f(u(x))$$

# $$\Longrightarrow \frac{1}{\Delta x} \left(h\left(x_{j + \frac{1}{2}}\right) -h\left(x_{j - \frac{1}{2}}\right)\right) = \left[f(u)\right]_x$$

# $$\Longrightarrow \frac{1}{\Delta x} \left(h\left(x_{j + \frac{1}{2}}\right) -h\left(x_{j - \frac{1}{2}}\right)\right) \approx \frac{1}{\Delta x}\left(\widehat{f}_{j + \frac{1}{2}} - \widehat{f}_{j - \frac{1}{2}}\right)$$

# $$\widehat{f}_{j + \frac{1}{2}} = h\left(x_{j + \frac{1}{2}}\right) = h_{j + \frac{1}{2}}$$

# $$\overline{h}_j^n = \frac{1}{\Delta x} \int_{x - \Delta x / 2}^{x + \Delta x / 2} h(s) \, ds = f\left(u_j^n\right)$$

# $$\left\{\overline{h}_j^n\right\}_j \overset{?}{\longrightarrow} \left\{h_{j + \frac{1}{2}}^n\right\}_j$$

# In[2]:


get_ipython().run_line_magic('matplotlib', 'inline')

import seaborn
import sympy

import weno_computations

seaborn.set_palette('husl')


# In[3]:


weno_computations.make_intro_plots(0)


# In[4]:


weno_computations.make_intro_plots(1)


# In[5]:


weno_computations.make_intro_plots(2)


# In[6]:


weno_computations.make_intro_plots(3)


# In[7]:


weno_computations.make_intro_plots(4)


# In[8]:


weno_computations.make_intro_plots(5)


# In[9]:


weno_computations.make_intro_plots(6)


# In[10]:


weno_computations.make_intro_plots(7)


# In[11]:


weno_computations.make_intro_plots(8)


# In[12]:


weno_computations.make_intro_plots(9)


# In[13]:


weno_computations.make_intro_plots(10)


# In[14]:


weno_computations.make_intro_plots(11)


# In[15]:


weno_computations.make_intro_plots(12)


# In[16]:


# The slide below can be produced with the following:
(approx_minus2, approx_minus1, 
 approx_zero, approx_all) = weno_computations.interp_simple_stencils()

approx_minus2 = approx_minus2.replace('=', '&=')
approx_minus1 = approx_minus1.replace('=', '&=')
approx_zero = approx_zero.replace('=', '&=')
approx_all = approx_all.replace('=', '&=')

all_math = r"""\begin{align*}
%s \\ %s \\ %s \\ %s \end{align*}
""" % (approx_minus2, approx_minus1, approx_zero, approx_all)


# $$u^{{(1)}}_{{j + \frac{1}{2}}} = \frac{1}{3}\overline{u}_{j-2} - \frac{7}{6}\overline{u}_{j-1} + \frac{11}{6}\overline{u}_{j}$$

# $$u^{{(2)}}_{{j + \frac{1}{2}}} = -\frac{1}{6}\overline{u}_{j-1} + \frac{5}{6}\overline{u}_{j} + \frac{1}{3}\overline{u}_{j+1}$$

# $$u^{{(3)}}_{{j + \frac{1}{2}}} = \frac{1}{3}\overline{u}_{j} + \frac{5}{6}\overline{u}_{j+1} - \frac{1}{6}\overline{u}_{j+2}$$

# $$u^{{(\ast)}}_{{j + \frac{1}{2}}} = \frac{1}{30}\overline{u}_{j-2} - \frac{13}{60}\overline{u}_{j-1} + \frac{47}{60}\overline{u}_{j} + \frac{9}{20}\overline{u}_{j+1} -\frac{1}{20}\overline{u}_{j+2}$$

# $$\frac{1}{3}\overline{u}_{j-2} - \frac{7}{6}\overline{u}_{j-1} + \frac{11}{6}\overline{u}_{j} = u\left(x_{j + \frac{1}{2}}\right) - \frac{\Delta x^3 u_0'''}{4} + \mathcal{O}\left(\Delta x^4\right)$$

# $$-\frac{1}{6}\overline{u}_{j-1} + \frac{5}{6}\overline{u}_{j} + \frac{1}{3}\overline{u}_{j+1} = u\left(x_{j + \frac{1}{2}}\right) + \frac{\Delta x^3 u_0'''}{12} + \mathcal{O}\left(\Delta x^4\right)$$

# $$\frac{1}{3}\overline{u}_{j} + \frac{5}{6}\overline{u}_{j+1} - \frac{1}{6}\overline{u}_{j+2} = u\left(x_{j + \frac{1}{2}}\right) - \frac{\Delta x^3 u_0'''}{12} + \mathcal{O}\left(\Delta x^4\right)$$

# $$u^{{(\ast)}}_{{j + \frac{1}{2}}} = u\left(x_{j + \frac{1}{2}}\right) - \frac{\Delta x^5 u_0'''}{60} + \mathcal{O}\left(\Delta x^6\right)$$

# $$\left(\Delta x\right) \overline{u} = \int_a^b u(x) \, dx = u_0(b - a) + u_0' \left(\frac{b^2 - a^2}{2}\right) + u_0'' \left(\frac{b^3 - a^3}{6}\right) + \cdots$$

# Putting the **W** in **WENO**:
# $$u^{{(\ast)}}_{{j + \frac{1}{2}}} = \frac{1}{10} u^{{(1)}}_{{j + \frac{1}{2}}} + \frac{6}{10} u^{{(2)}}_{{j + \frac{1}{2}}} + \frac{3}{10} u^{{(3)}}_{{j + \frac{1}{2}}}$$

# What about the other letters?

# What about shocks?

# In[17]:


weno_computations.discontinuity_to_volume_single_cell(0)


# In[18]:


weno_computations.discontinuity_to_volume_single_cell(1)


# In[19]:


weno_computations.make_shock_plot_single_cell()


# $$\frac{1}{10} u^{{(1)}}_{{j + \frac{1}{2}}} + \frac{6}{10} u^{{(2)}}_{{j + \frac{1}{2}}} + \frac{3}{10} u^{{(3)}}_{{j + \frac{1}{2}}} \longrightarrow \frac{2}{3} u^{{(2)}}_{{j + \frac{1}{2}}} + \frac{1}{3} u^{{(3)}}_{{j + \frac{1}{2}}}$$

# In[20]:


weno_computations.discontinuity_to_volume()


# In[21]:


weno_computations.make_shock_plot()


# Smoothness indicator when approximating a point in the interval $I = \left(x_j - \frac{\Delta x}{2}, x_j + \frac{\Delta x}{2}\right]$

# $$\beta = \sum_k \Delta x^{2k - 1} \int_I \left(\frac{d^k}{dx^k} p\right)^2 dx$$

# Large $\beta$ indicates that the approximation $p(x)$ is not smooth nearby our point.

# Recall the **W** in **WENO**:
# $$u^{{(\ast)}}_{{j + \frac{1}{2}}} = \frac{1}{10} u^{{(1)}}_{{j + \frac{1}{2}}} + \frac{6}{10} u^{{(2)}}_{{j + \frac{1}{2}}} + \frac{3}{10} u^{{(3)}}_{{j + \frac{1}{2}}}$$

# Rather than using these directly, incorporate a penalty for each polynomial approximation's (lack of) smoothness: $\beta^{(1)}$, $\beta^{(2)}$, $\beta^{(3)}$.

# $$\widetilde{w}^{(1)} = \frac{0.1}{\left(\varepsilon + \beta^{(1)}\right)^2}, \;
# \widetilde{w}^{(2)} = \frac{0.6}{\left(\varepsilon + \beta^{(2)}\right)^2}, \;
# \widetilde{w}^{(3)} = \frac{0.3}{\left(\varepsilon + \beta^{(3)}\right)^2}$$

# $$w^{(j)} = \frac{\widetilde{w}^{(j)}}{\widetilde{w}^{(1)} + \widetilde{w}^{(2)} + \widetilde{w}^{(3)}}$$

# $$w^{(1)} u^{{(1)}}_{{j + \frac{1}{2}}} + w^{(2)} u^{{(2)}}_{{j + \frac{1}{2}}} + w^{(3)} u^{{(3)}}_{{j + \frac{1}{2}}}$$

# <h1 class="bigtitle">Back to Conservation Law</h1>

# $$\dot{\overline{u}} = -\frac{1}{\Delta x}\left(\widehat{f}\left(u_{j + \frac{1}{2}}^-, u_{j + \frac{1}{2}}^+\right) - \widehat{f}\left(u_{j - \frac{1}{2}}^-, u_{j - \frac{1}{2}}^+\right)\right)$$

# $$u_{j + \frac{1}{2}}^- \longrightarrow \text{Reconstruct from } \overline{u}_{j - 2}, \overline{u}_{j - 1}, \overline{u}_{j}, \overline{u}_{j + 1}, \overline{u}_{j + 2}$$

# $$u_{j + \frac{1}{2}}^+ \longrightarrow \text{Reconstruct from } \overline{u}_{j - 1}, \overline{u}_{j}, \overline{u}_{j + 1}, \overline{u}_{j + 2}, \overline{u}_{j + 3}$$

# $$\dot{\overline{u}} = L\left(\overline{u}\right)$$

# TV-diminishing / strong stability preserving RK scheme
# 
# $$\begin{array}{c | c c c}
#   0 &     &     &   \\
#   1 &   1 &     &   \\
# 1/2 & 1/4 & 1/4 &   \\
# \hline
#     & 1/6 & 1/6 & 2/3
# \end{array}$$

# ![Stability Preserving](tvd_scheme.png "TVD RK Scheme")

# <h1 class="bigtitle">Example</h1>

# ![WENO Example Part 1](traffic_flow_example1.png "WENO Traffic Flow")

# ![WENO Example Part 2](traffic_flow_example2.png "WENO Traffic Flow")