Make Movies as Numerical Solutions of Partial Differential Equations

[Eli]

This is a snapshot of numerical solutions of Partial Differential Equations - PDEs. We particularly want to show how to numerically solve PDEs of the form \(u_{t} + f(u)_{x} = 0\), where \( u = u(x,t)\) is a function of space, \(x\) and time, \(t\). The PDEs of this form usually describe transport phenomena, and are referred to as advection equations. Categorically, hyperbolic conservation laws (HCLs) fall into this group. Non-linear HCLs are in most cases difficult to solve. One major problem is that, HCLs have a characteristic of developing discontinuities or shock waves at a later time, as they evolve, even if their initial conditions are smooth.

Simple numerical schemes such as Upwind can be used to solve such differential equations, but optimal solutions can be achieved by combining various schemes. You can consult the post High-order Relaxation Approaches for Adjoint-based Optimization to access exhaustive work on HCLs and their solutions using a combination of a number of approaches.

We will consider inviscid Burgers' equation (see full details here) for this demo.

Requirements

You need a working installation of python-scitools. Python Scitools works best with Python 2.7.x.

If you are using Debian, or a Debian based platform like Ubuntu-Linux, you can install scitools using sudo:

$sudo apt-get install python-scitools

Be cautious, that if you have another version of Python installed (say Python 3.6.x) and use it as a default Python, this way of installing Scitools may not work because it may create conflicts between the two versions, or likely to get a message "No module named scitools" when you run a code with scitools import. However, Scitools comes with a lot of other packages, system wide installation is not recommended since this may break your system and render parts of it unusable due to changes that Scitools will make.

You can make scitools install by specifying location:

$python setup.py install --prefix=/path/to/install/directory

Then setup the PYTHONPATH to include the path where the package is installed, in the ~/.bashrc (.bash_profile), like,

export PYTHONPATH=/path/to/install/directory:$PYTHONPATH

This will enable Python and other packages to find scitools from wherever folder you have installed it.

If you only have Python 2.7 installed in your system using Anaconda, you can use conda package to install scitools:

$conda install --channel johannr scitools

The best ever and safe way to install scitools is to use virtual environments - virtualenv.

Using Virtual Python Environments

Create a Python Virtual environment, let's assume you name it "PDEs" and install scitools and any other required packages, mostly numpy.

Create folder for virtualenv:

$mkdir Enviroments

Move inside Environments folder:

$cd Environments

Make virtualenv named PDEs:

$virtualenv PDEs

Activate the virtualenv:

$source PDEs/bin/activate

Install requirements - numpy and scitools:

$pip install numpy
$git clone https://github.com/hplgit/scitools.git
$cd scitools
$python setup.py install

If your system uses Python 2.7, you do not need to re-install it inside virtualenv, instead create a virtualenv PDEs on fly by specifying the python version to be used with it. In our case, we want to use python 2.7, so we would run the following command:

$virtualenv -p /usr/bin/python2.7 PDEs

This is especially very important and useful if you have several Python versions installed in your system, to avoid possible problems.

Next, make directory where you will put your code, in our case, we will make directory a directory under /home/user/ called PDEs_Simulations:

$mkdir PDEs_Simulations

Create a Python file, for example, named "burgers.py":

$touch burgers.py

and then open it with your favourite editor/IDE to develop your code in it or cut/copy and paste. Make sure you run the code while inside the virtual environment PDEs.

You can choose between Python Matplotlib, Gnuplot, Pylab, etc., as plotting backends, see customization options.

By default, plotting is done by Easyviz, scitools.easyviz module - a federated interface to various packages for scientific plotting and visualization.

Terms of discrete PDEs, such as \(u_{k-1}, \ u_{k}, \ u_{k+1}\) are normally, respectively, translated to \(u[k-1], \ u[k], \ u[k+1]\) in numerical Python code.

You can then import scitools as

$from scitools.std import *

or minimally in this way,

import scitools.std as st

and specify it before any plotting command: st.command.

For sure, there are lot more we can do with Python Scitools such as decorating and adding controls to our plots, making movies and animations in different formats. Ability to make movies in different formats depends on the decoders installed in your system.

Below is a minimal working example implementation of the Burgers' equation as introduced in this topic:

Don't forget to clean a lot of image files that will be created successively using:

Code: [Select all] [Expand/Collapse]

1. for filename in glob.glob('/home/user/PDEs_Simulations/wave*.png'):
2.     os.remove(os.path.join("/home/user/PDEs_Simulations", filename))

The path /home/user/PDEs_Simulations depends on the location of your PDEs_Simulations directory.

Code: [Select all] [Expand/Collapse]

1. """Solve Burger's Equation and
2. simulate its evolution with time"""
3. import numpy as np
4. from scitools.std import *
5. import time
6. import glob, os
7.  
8. N = 850 #Simulate for long enough time
9. M = 400
10. dx = 2 * np.pi / M
11. cfl = 0.55
12. a = 1.3
13. dt = cfl * dx / np.sqrt(a)
14.  
15. #Initial conditions
16. def initialize():
17.     return 0.5 + np.sin(x)
18.  
19. x = np.linspace(0.0, 2 * np.pi, M + 1)
20. u = initialize()
21.  
22. def evolve_in_One_Timestep(u):
23.     #Create an array to store the new values of the solution
24.     u_new = np.zeros((M + 1,), float)
25.     for k in range(1, M):  # M-1
26.         u_new[k] = u[k] - 0.25 * (u[k + 1] ** 2 - u[k - 1] ** 2) * (dt / dx) \
27.         + 0.5 * np.sqrt(a) * (u[k + 1] - 2 * u[k] + u[k - 1]) * (dt / dx)
28.     #Apply boundary conditions
29.     u_new[0] = u_new[1]
30.     u_new[M] = u_new[M - 1]
31.     return u_new
32.  
33. iteration = 0
34. evolve = initialize()
35. for n in range(0, N):
36.     plot(np.linspace(0.0, 1.0, M + 1), evolve, 'r', linewidth=1.0, show=False,
37.          xlabel='x', ylabel='u(x,t)', title="Evolution of wave Equation", \
38.          markersize=8, legend='step=%4.0f' % n,
39.          hardcopy="wave%04d.png" % iteration)
40.     evolve = evolve_in_One_Timestep(evolve)
41.     iteration += 1
42.  
43. #Make movie
44. movie("wave*.png")
45.  
46. #Explicitly specify .gif movie using controls, fps means # of frames per second
47. movie('wave*.png', encoder='convert', fps=12,
48.       output_file='Burgers.gif')
49.  
50. #Generate an avi movie
51. movie('wave*.png', encoder='ffmpeg', fps=12,
52.       output_file='Burgers.avi')
53.  
54. savefig("Shock_wave.png")
55.  
56. #Clean all created images
57. for filename in glob.glob('/home/tssfl/PDEs_Simulations/wave*.png'):
58.     os.remove(os.path.join("/home/tssfl/PDEs_Simulations", filename))

Outputs should be the files movie.gif, Burgers.gif, Burgers.avi and Shock_wave.png.

[Eli]

You can now git clone the code from GitHub and run it right away:

https://github.com/elimboto/pde_simulations