<jupyter_start><jupyter_text>Text provided under Creative Commons Attribution License CC-BY.  All Code is made availableunder the FSF-approved BSD-3 license (c) <NAME>, <NAME> 2017. Thanks toNSF for support via CAREER award 1149784.   CFDPython_4.pySource: https://github.com/barbagroup/CFDPython/blob/master/lessons/03_CFL_Condition.ipynbThis if from a class on learning to program CFD using Python by Dr. <NAME> Step 4: 1-D Burgers' Equation - convection and diffusionAlso called viscous Burgers' equation for convection and diffusion in one dimension.  You can read about Burgers' Equation on its wikipedia page.The equation looks like this:As the equation shows, it is a combination of non-linear convection and diffusion equations seen separately in the previous lessons 1-3.  We can solve this equation by using a taylor series as in the previous lessons.  Using forward difference for time, backward difference for space and our 2nd-order method for the second derivitive yields:We can solve for our only unknow , and step in time as follows:Initial and Boundary ConditionsTo investigate some interesting properties of Berger's equation, it is helpful to use different initial and boundary conditions than we used in our previous steps.  Our initial conditions are going to be<img class="math math-display" alt="$$\begin{eqnarray}u &amp;amp;=&amp;amp; -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\\phi &amp;amp;=&amp;amp; \exp \bigg(\frac{-x^2}{4 \nu} \bigg) + \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu} \bigg)\end{eqnarray}$$" src="https://render.githubusercontent.com/render/math?math=%5Cbegin%7Beqnarray%7D%0Au%20%26amp%3B%3D%26amp%3B%20-%5Cfrac%7B2%20%5Cnu%7D%7B%5Cphi%7D%20%5Cfrac%7B%5Cpartial%20%5Cphi%7D%7B%5Cpartial%20x%7D%20%2B%204%20%5C%5C%5C%0A%5Cphi%20%26amp%3B%3D%26amp%3B%20%5Cexp%20%5Cbigg%28%5Cfrac%7B-x%5E2%7D%7B4%20%5Cnu%7D%20%5Cbigg%29%20%2B%20%5Cexp%20%5Cbigg%28%5Cfrac%7B-%28x-2%20%5Cpi%29%5E2%7D%7B4%20%5Cnu%7D%20%5Cbigg%29%0A%5Cend%7Beqnarray%7D&amp;mode=display">This has an analytical solution given by<img class="math math-display" alt="$$\begin{eqnarray}u &amp;amp;=&amp;amp; -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\\phi &amp;amp;=&amp;amp; \exp \bigg(\frac{-(x-4t)^2}{4 \nu (t+1)} \bigg) + \exp \bigg(\frac{-(x-4t -2 \pi)^2}{4 \nu(t+1)} \bigg)\end{eqnarray}$$" src="https://render.githubusercontent.com/render/math?math=%5Cbegin%7Beqnarray%7D%0Au%20%26amp%3B%3D%26amp%3B%20-%5Cfrac%7B2%20%5Cnu%7D%7B%5Cphi%7D%20%5Cfrac%7B%5Cpartial%20%5Cphi%7D%7B%5Cpartial%20x%7D%20%2B%204%20%5C%5C%5C%0A%5Cphi%20%26amp%3B%3D%26amp%3B%20%5Cexp%20%5Cbigg%28%5Cfrac%7B-%28x-4t%29%5E2%7D%7B4%20%5Cnu%20%28t%2B1%29%7D%20%5Cbigg%29%20%2B%20%5Cexp%20%5Cbigg%28%5Cfrac%7B-%28x-4t%20-2%20%5Cpi%29%5E2%7D%7B4%20%5Cnu%28t%2B1%29%7D%20%5Cbigg%29%0A%5Cend%7Beqnarray%7D&amp;mode=display"> Our boundary condition will be:This is called a periodic boundary condition.Saving Time with SymPyThe initial condition we're using for Burgers' Equation can be a bit of a pain to evaluate by hand. The derivative $\frac{\partial \phi}{\partial x}$ isn't too terribly difficult, but it would be easy to drop a sign or forget a factor of $x$ somewhere, so we're going to use SymPy to help us out.SymPy is the symbolic math library for Python. It has a lot of the same symbolic math functionality as Mathematica with the added benefit that we can easily translate its results back into our Python calculations (it is also free and open source).Start by loading the SymPy library, together with our favorite library, NumPy.  Then we're going to tell SymPy that we want all of its output to be rendered using $\LaTeX$. This will make our Notebook beautiful!<jupyter_code>#Remember: Comments in python are denoted by the pound sign 
import numpy                   # Load the numpy maxrix algebra library
import sympy                   # To simplify our mathematics use this symbolic math library.

#tell Sympy that we want our output rendered using LATEX.  This will make our notebook beautiful?!
from sympy import init_printing
init_printing(use_latex=True)<jupyter_output><empty_output><jupyter_text>Now we set up the symbolic variable we will use in sympy to compute phi and output the equation for phi in LATEX format to verify we got what we thought we were programming.<jupyter_code>x, nu, t = sympy.symbols('x nu t')
phi = (sympy.exp(-(x - 4 * t)**2 / (4 * nu * (t + 1))) +
       sympy.exp(-(x - 4 * t - 2 * sympy.pi)**2 / (4 * nu * (t + 1))))
phi<jupyter_output><empty_output><jupyter_text>The above output is small and hard to read, but corresponds to our equation for phi.  Nowset up our equation for the partial derivative of phi with respect to X.<jupyter_code>phiprime = phi.diff(x)
phiprime<jupyter_output><empty_output><jupyter_text>To see the unrendered version use the Python print command<jupyter_code>print(phiprime)<jupyter_output>-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))<jupyter_text>Now What?Now we have a Python version of phi and it's derivative.  Next we want to write out the full initial condition equation and then translate it into a usable python expression.  For this we will use the lambdify function from the SymPy.utilities library which takes a SymPy symbolic equation and turns it into a callable function."Numerical" VS "Analytical" solutions"Numerical" usually indicates an approximate solution obtained by methods of numerical analysis.  "Analytical" solutions are exact and obtained by methods of symbolic manipulation, derived using analysis. Typically Analytical solutions are signifitantly more compute intensive than Numerical methods, but provide exact solutions as opposed to the numerical solutions that generally involve some approximation to simplify their calculations.  Below we are going to compare a Numberical and Analytical solution to see what the Numeric approximation error looks like.<jupyter_code>from sympy.utilities.lambdify import lambdify

#initial condition equation:
u = -2 * nu *(phiprime / phi) + 4
print (u)<jupyter_output>-2*nu*(-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)))/(exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1))) + exp(-(-4*t + x)**2/(4*nu*(t + 1)))) + 4<jupyter_text>Lambdify:The sympy.utilities.lambdify module provides convenient functions to transform sympy expressions to lambda functions which can be used to calculate numerical values very fast. In addition, it converts SymPy names to the names used in the given numerical library, usually NumPy. To lambdify an expression into a usable function, we tell lambdify which variable to request and the function we want to plug them into.<jupyter_code>ufunc = lambdify((t, x, nu), u)
print(ufunc(1,4,3))<jupyter_output>3.49170664206445<jupyter_text>Back to Burgers' EquationNow that we have the initial conditions set, we can finish setting up the problem. Once the setup is complete, we will plot the initial conditions using our lambdify-ed function.<jupyter_code>#load our plotting routine and tell it to output our plots inline with our
# Junyper Notebook, rather than in a separate window.
from matplotlib import pyplot
%matplotlib inline

#Now declare our variables
nx = 101
nt = 100
dx = 2 * numpy.pi / (nx-1)
nu = 0.07
dt = dx * nu

x = numpy.linspace(0, 2 * numpy.pi, nx)
un = numpy.empty(nx)
t = 0

u = numpy.asarray([ufunc(t, x0, nu) for x0 in x])
#Output the resulting numerical array.
u
#plot the outline of our figure.
pyplot.figure(figsize=(11,7), dpi=100)
#set the x-axis to have markers
pyplot.plot(x, u, marker='o', lw=2)
#set the x-axis limits from 0 to 2*pi
pyplot.xlim({0, 2 * numpy.pi})
#set the y-axis limits
pyplot.ylim([0,10]);<jupyter_output><empty_output><jupyter_text>This is a numeric solution to Burgers' equation.  This is definitely not the hat function we were previously dealing with until now.  This function is called a "saw-tooth function".  Lets proceed and see what happens.Periodic Boundary ConditionsOne of the big differences between Step 3 and the previous lessons is the use of periodic boundary conditions. If you experimented with step 1 and 2 and make the simulation run long enough (by increasing nt) you will notice the wave will keep moving to the right until it no longer even shows up on the plot.  With periodic boundary conditions, when a point gets to the right-hand side of the frame, it wraps around back to the front of the frame.  Recall the discretization that we worked out at the beginning of this notebook:"What does"mean " when "" is already at the end of the frame?"Think about this for a minute before proceeding.Below we are going to compute the analytical solution and plot the solutions (the Numerical or Computational and Analytical) together.<jupyter_code>for n in range(nt):
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] - un[i] * dt / dx *(un[i] - un[i-1]) + nu * dt / dx**2 *\
                (un[i+1] - 2 * un[i] + un[i-1])
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu * dt / dx**2 *\
                (un[1] - 2 * un[0] + un[-2])
    u[-1] = u[0]
        
u_analytical = numpy.asarray([ufunc(nt * dt, xi, nu) for xi in x])
pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.plot(x,u, marker='o', lw=2, label='Computational')
pyplot.plot(x, u_analytical, label='Analytical')
pyplot.xlim([0, 2 * numpy.pi])
pyplot.ylim([0, 10])
pyplot.legend();<jupyter_output><empty_output><jupyter_text>What Next?The subsequent steps from 4 to 11 will extend the 1D finite-difference formulas (diffusion, convection and burgers' equations) into two dimensions.  Just apply the definition -- a partial derivative with respect to x is the variation in the x direction while keeping y constant.Before moving on, take a slight break and learn about array operations with mumpy.<jupyter_code>from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11817102137767221, 'fraction_numerical': 0.0447347585114806, 'mean_word_length': 4.8955659276546095, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '660050', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3603, 'n_tokens_neox': 3356, 'n_words': 1195, 'path': "lessons/05_Step_4-1D_Burger's_Equation-Convection_and_Diffusion.ipynb", 'repo_name': 'BBFlyer1/CFDPython', 'size': '10134'}