Merge branch 'master' into dse_fix

.travis.yml

@@ -55,5 +55,5 @@ script:
   - python examples/benchmark.py test -P acoustic -a
   - python examples/seismic/acoustic/acoustic_example.py
   - py.test -vs examples/seismic/tutorials
-  - py.test -vs examples/diffusion/example_diffusion.py
+  - py.test -vs examples/cfd
   - if [[ $DEVITO_ARCH == 'gcc-5' ]]; then ./docs/deploy.sh; fi

docs/deploy.sh

@@ -45,9 +45,10 @@ git config user.name "Travis CI"
 git config user.email "$COMMIT_AUTHOR_EMAIL"
 
 # If there are no changes to the compiled out (e.g. this is a README update) then just bail.
-if git diff --quiet; then
-    echo "No changes to the output on this push; exiting."
-    exit 0
+# NOTE: Disabled since it also suppresses re-compiled tutorial links, etc.
+#if git diff --quiet; then
+#    echo "No changes to the output on this push; exiting."
+#    exit 0
 fi
 
 # Commit the "changes", i.e. the new version.

docs/tutorials.rst

@@ -12,14 +12,34 @@ by checking out our the repository on github_
 Computational Fluid Dynamics
 ----------------------------
 
-**COMING SOON:** A tutorial series that introduces the core features
-of Devito through a set of classic examples from Computational Fluid
-Dynamics (CFD). The series is based on the excellent tutorial blog
-`CFD Python: 12 steps to Navier-Stokes
+A tutorial series that introduces the core features of Devito through
+a set of classic examples from Computational Fluid Dynamics (CFD). The
+series is based on the excellent tutorial blog `CFD Python: 12 steps
+to Navier-Stokes
 <http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/>`_
 by the Lorena A. Barba Group and focusses on the implementation with
 Devito rather than pure CFD or finite difference theory.
 
+* `01 - Linear convection
+  <http://nbviewer.jupyter.org/github/opesci/devito/blob/master/examples/cfd/test_01_convection.ipynb>`_
+  - Building a linear operator with simple zero BCs.
+* `02 - Nonlinear convection
+  <http://nbviewer.jupyter.org/github/opesci/devito/blob/master/examples/cfd/test_02_convection_nonlinear.ipynb>`_
+  - Building an operator with coupled equations and simple BCs.
+* `03 - Diffusion
+  <http://nbviewer.jupyter.org/github/opesci/devito/blob/master/examples/cfd/test_03_diffusion.ipynb>`_
+  - Building an second-order operator with simple BCs.
+* `04 - Burgers' equation
+  <http://nbviewer.jupyter.org/github/opesci/devito/blob/master/examples/cfd/test_04_burgers.ipynb>`_
+  - Coupled operator with mixed discretizations and simple BCs.
+* `05 - Laplace equation
+  <http://nbviewer.jupyter.org/github/opesci/devito/blob/master/examples/cfd/test_05_laplace.ipynb>`_
+  - Steady-state example with Python convergence loop and Neumann BCs.
+* `06 - Poisson equation
+  <http://nbviewer.jupyter.org/github/opesci/devito/blob/master/examples/cfd/test_06_poisson.ipynb>`_
+  - Pseudo-timestepping example with kernel-driven diffusion loop.
+
+
 Seismic Modelling and Inversion
 -------------------------------
 

examples/cfd/__init__.py

@@ -0,0 +1 @@
+from .tools import *  # noqa

examples/diffusion/example_diffusion.py → examples/cfd/example_diffusion.py

@@ -1,7 +1,16 @@
 """
+NOTE: This is a legacy example. For more up-to-date examples
+please see the tutorials under `examples/cfd/`.
+
 This example encodes multiple ways to solve the 2D diffusion equations
 using an explicit finite difference scheme with fixed boundary values
-and a given initial value for the density.
+and a given initial value for the density. The example also includes
+demo implementations using more traditional methods, including pure
+Python, vectorized NumPy and a vectorized symbolic implementation
+using SymPy's `lambdify()` functionality. For a visual demonstration
+and comparison between pure Python and NumPy, for example, run:
+
+python example_diffusion.py run -m python numpy --show
 """
 import time
 from argparse import ArgumentParser

examples/cfd/tools.py

@@ -0,0 +1,50 @@
+from mpl_toolkits.mplot3d import Axes3D  # noqa
+
+import numpy as np
+from matplotlib import pyplot, cm
+
+
+def plot_field(field, xmax=2., ymax=2., zmax=None, view=None, linewidth=0):
+    """Utility plotting routine for 2D data
+
+    :param field: Numpy array with field data to plot
+    :param xmax: (Optional) Length of the x-axis
+    :param ymax: (Optional) Length of the y-axis
+    :param view: (Optional) View point to intialise
+    """
+    x_coord = np.linspace(0, xmax, field.shape[0])
+    y_coord = np.linspace(0, ymax, field.shape[1])
+    fig = pyplot.figure(figsize=(11, 7), dpi=100)
+    ax = fig.gca(projection='3d')
+    X, Y = np.meshgrid(x_coord, y_coord)
+    ax.plot_surface(X, Y, field[:], cmap=cm.viridis, rstride=1, cstride=1,
+                    linewidth=linewidth, antialiased=False)
+
+    # Enforce axis measures and set view if given
+    ax.set_xlim(0., xmax)
+    ax.set_ylim(0., ymax)
+    if zmax is not None:
+        ax.set_zlim(1., zmax)
+    if view is not None:
+        ax.view_init(*view)
+
+    # Label axis
+    ax.set_xlabel('$x$')
+    ax.set_ylabel('$y$')
+
+    pyplot.show()
+
+
+def init_hat(field, dx, dy, value=2., bgvalue=1.):
+    """Set "hat function" initial condition on an array:
+
+    u(.5<=x<=1 && .5<=y<=1 ) is 2
+
+    :param field: Numpy array with field data to plot
+    :param dx: Spacing in the x-dimension
+    :param dy: Spacing in the y-dimension
+    :param value: Value of the top part of the function, default=2.
+    :param bgvalue: Background value for the bottom of the function, default=1.
+    """
+    field[:] = bgvalue
+    field[int(.5 / dx):int(1 / dx + 1), int(.5 / dy):int(1 / dy + 1)] = value

examples/seismic/acoustic/operators.py

@@ -49,7 +49,7 @@ def ForwardOperator(model, source, receiver, time_order=2, space_order=4,
     #   u[ti + 1, ...] - accesses the forward stencil value
     ti = u.indices[0]
     src_term = src.inject(field=u, u_t=ti + 1, offset=model.nbpml,
-                          expr=src * dt * dt / m, p_t=time)
+                          expr=src * dt**2 / m, p_t=time)
 
     # Create interpolation expression for receivers
     rec_term = rec.interpolate(expr=u, u_t=ti, offset=model.nbpml)
@@ -95,7 +95,7 @@ def AdjointOperator(model, source, receiver, time_order=2, space_order=4, **kwar
     # Construct expression to inject receiver values
     ti = v.indices[0]
     receivers = rec.inject(field=v, u_t=ti - 1, offset=model.nbpml,
-                           expr=rec * dt * dt / m, p_t=time)
+                           expr=rec * dt**2 / m, p_t=time)
 
     # Create interpolation expression for the adjoint-source
     source_a = srca.interpolate(expr=v, u_t=ti, offset=model.nbpml)