Merge pull request #100 from devitocodes/devito

Contributing guidelines, corrections and tests.

Author: ggorman

CONTRIBUTING.md

@@ -775,7 +775,7 @@ u^{n+1}_i = u^{n-1}**i - C(u**{i+1}^n-u_{i-1}^n)\tp
 $$
 A special scheme is needed to compute $u^1$, but we leave that problem for
 now. Anyway, this special scheme can be found in
-[`advec1D.py`](https://github.com/hplgit/fdm-book/tree/master/src/advec/advec1D.py).
+[`advec1D.py`](https://github.com/devitocodes/devito_book/tree/main/src/advec/advec1D.py).
 
 ### Implementation
 We now need to work with three time levels and must modify our solver a bit:
@@ -1171,7 +1171,7 @@ def run(scheme='UP', case='gaussian', C=1, dt=0.01):
     print 'Integral of u:', integral.max(), integral.min()
 ```
 The complete code is found in the file
-[`advec1D.py`](https://github.com/hplgit/fdm-book/tree/master/src/advec/advec1D.py).
+[`advec1D.py`](https://github.com/devitocodes/devito_book/tree/main/src/advec/advec1D.py).
 
 ## A Crank-Nicolson discretization in time and centered differences in space {#sec-advec-1D-CN}
 

chapters/advec/advec.qmd

@@ -226,7 +226,7 @@ $$ {#eq-form-Dop-tn3-fw}
 
 Application of finite difference operators to polynomials and exponential
 functions, resulting in the formulas above, can easily be computed by
-some `sympy` code (from the file [`lib.py`](https://github.com/hplgit/fdm-book/tree/master/src/formulas/lib.py)):
+some `sympy` code (from the file `src/formulas/lib.py` in this repository):
 
 ```python
 from sympy import *

chapters/appendices/formulas/formulas.qmd

@@ -36,24 +36,19 @@ $$
 [D_{2x}u]^n_i=0,
 $$
 at a boundary point $i$. The details of how the numerical scheme
-is worked out are described in ref[Sections @sec-wave-pde2-Neumann
-and @sec-wave-pde2-var-c][ in [@Langtangen_deqbook_wave]][in
-the document [Finite difference methods for wave motion](http://tinyurl.com/k3sdbuv/pub/wave)
-[@Langtangen_deqbook_wave]].
+is worked out are described in @sec-wave-pde2-Neumann
+and @sec-wave-pde2-var-c.
 
 ## A solver function
 
 The general initial-boundary value problem
 solved by finite difference methods can be implemented as shown in
 the following `solver` function (taken from the
-file [`wave1D_dn_vc.py`](https://github.com/hplgit/fdm-book/tree/master/src/wave/wave1D/wave1D_dn_vc.py)).
+file [`wave1D_dn_vc.py`](https://github.com/devitocodes/devito_book/tree/main/src/wave/wave1D/wave1D_dn_vc.py)).
 This function builds on
 simpler versions described in
-ref[Sections @sec-wave-pde1-impl, @sec-wave-pde1-impl-vec
-@sec-wave-pde2-Neumann, and @sec-wave-pde2-var-c][ in
-[@Langtangen_deqbook_wave]][in
-the document [Finite difference methods for wave motion](http://tinyurl.com/k3sdbuv/pub/wave)
-[@Langtangen_deqbook_wave]].
+@sec-wave-pde1-impl, @sec-wave-pde1-impl-vec,
+@sec-wave-pde2-Neumann, and @sec-wave-pde2-var-c.
 There are several quite advanced
 constructs that will be commented upon later.
 The code is lengthy, but that is because we provide a lot of
@@ -286,7 +281,7 @@ Such actions must be
 taken care of outside the `solver` function, more precisely in the
 `user_action` function that is called at every time level.
 
-We have, in the [`wave1D_dn_vc.py`](https://github.com/hplgit/fdm-book/tree/master/src/wave/wave1D/wave1D_dn_vc.py)
+We have, in the [`wave1D_dn_vc.py`](https://github.com/devitocodes/devito_book/tree/main/src/wave/wave1D/wave1D_dn_vc.py)
 code, implemented the `user_action`
 callback function as a class `PlotAndStoreSolution` with a
 `__call__(self, x, t, t, n)` method for the `user_action` function.
@@ -357,7 +352,7 @@ The `hashed_input` argument, used to name the
 resulting archive file with all solutions, is supposed to be a
 hash reflecting all import parameters in the problem such that this
 simulation has a unique name.
-The `hashed_input` string is made in the [`solver`](https://github.com/hplgit/fdm-book/tree/master/src/wave/wave1D/wave1D_dn_vc.py) function, using the `hashlib`
+The `hashed_input` string is made in the [`solver`](https://github.com/devitocodes/devito_book/tree/main/src/wave/wave1D/wave1D_dn_vc.py) function, using the `hashlib`
 and `inspect` modules, based on the arguments to `solver`:
 
 ```python
@@ -419,10 +414,7 @@ solutions to machine precision. With Dirichlet boundary conditions
 we can construct a function that is linear in $t$ and quadratic in
 $x$ that is also an exact solution of the scheme, while with Neumann
 conditions we are left with testing just a constant solution
-(see comments in ref[Section @sec-wave-pde1-verify][ in
-[@Langtangen_deqbook_wave]][in
-the document [Finite difference methods for wave motion](http://tinyurl.com/k3sdbuv/pub/wave)
-[@Langtangen_deqbook_wave]]).
+(see comments in @sec-wave-pde1-verify).
 
 ### Convergence rates
 
@@ -435,9 +427,7 @@ $$
 r = \frac{\ln(E_{i}/E_{i-1})}{\ln(h_{i}/h_{i-1})},
 $$
 where $E_i$ is the error corresponding to $h_i$ and $E_{i-1}$ corresponds to
-$h_{i-1}$. ref[Section @sec-wave-pde2-fd-standing-waves][ in [@Langtangen_deqbook_wave]][The
-section "Using an analytical solution of physical significance": ""
-in [@Langtangen_deqbook_wave]] explains the details of this type of verification and how
+$h_{i-1}$. @sec-wave-pde2-fd-standing-waves explains the details of this type of verification and how
 we introduce the single discretization parameter $h=\Delta t = \hat c\Delta t$,
 for some constant $\hat c$. To compute the error, we had to rely on
 a global variable in the user action function. Below is an implementation
@@ -449,15 +439,12 @@ error (which is always considered an advantage).
 ```
 
 The returned sequence `r` should converge to 2 since the error
-analysis in ref[Section @sec-wave-pde1-analysis][in
-[@Langtangen_deqbook_wave]][the section "Analysis of the difference
-equations": "" in [@Langtangen_deqbook_wave]] predicts various error measures to behave
+analysis in @sec-wave-pde1-analysis predicts various error measures to behave
 like $\Oof{\Delta t^2} + \Oof{\Delta x^2}$. We can easily run
 the case with standing waves and the analytical solution
 $u(x,t) = \cos(\frac{2\pi}{L}t)\sin(\frac{2\pi}{L}x)$. The call will
 be very similar to the one provided in the `test_convrate_sincos` function
-in ref[Section @sec-wave-pde1-impl-verify-rate][ in [@Langtangen_deqbook_wave]][the section "Verification: convergence rates": "" in
-[@Langtangen_deqbook_wave]], see the file [`wave1D_dn_vc.py`](https://github.com/hplgit/fdm-book/tree/master/src/wave/wave1D/wave1D_dn_vc.py) for details.
+in @sec-wave-pde1-impl-verify-rate, see the file `src/wave/wave1D/wave1D_dn_vc.py` for details.
 
 Many who know about class programming prefer to organize their software
 in terms of classes. This gives a richer application programming interface
@@ -482,7 +469,7 @@ how this may be counteracted by introducing a super class `Parameters` that allo
 code to be parameterized. In addition, it is convenient to collect the
 arrays that describe the mesh in a special `Mesh` class and make
 a class `Function` for a mesh function (mesh point values and its mesh).
-All the following code is found in [`wave1D_oo.py`](https://github.com/hplgit/fdm-book/tree/master/src/softeng2/wave1D_oo.py).
+All the following code is found in [`wave1D_oo.py`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave1D_oo.py).
 
 ## Class Parameters
 
@@ -1222,22 +1209,21 @@ is reproduced (within machine precision).
 
 ## Speeding up Cython code {#sec-wave2D3D-impl-Cython}
 
-We now consider the [`wave2D_u0.py`](https://github.com/hplgit/fdm-book/tree/master/src/wave/wave2D_u0/wave2D_u0.py)
+We now consider the [`wave2D_u0.py`](https://github.com/devitocodes/devito_book/tree/main/src/wave/wave2D_u0/wave2D_u0.py)
 code for solving the 2D linear wave equation with constant wave
 velocity and homogeneous Dirichlet boundary conditions $u=0$.
 We shall in the present chapter extend this code with computational
 modules written in other languages than Python. This extended version is
-called [`wave2D_u0_adv.py`](https://github.com/hplgit/fdm-book/tree/master/src/softeng2/wave2D_u0_adv.py).
+called [`wave2D_u0_adv.py`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave2D_u0_adv.py).
 
 The `wave2D_u0.py` file contains a `solver` function, which calls an
 `advance_*` function to advance the numerical scheme one level forward
 in time.  The function `advance_scalar` applies standard Python loops
 to implement the scheme, while `advance_vectorized` performs
 corresponding vectorized arithmetics with array slices. The statements
-of this solver are explained in ref[Section @sec-wave-2D3D-impl, in
-particular Sections @sec-wave2D3D-impl-scalar and
-@sec-wave2D3D-impl-vectorized][ in [@Langtangen_deqbook_wave]][in
-the document [Finite difference methods for wave motion](http://tinyurl.com/k3sdbuv/pub/wave) [@Langtangen_deqbook_wave]].
+of this solver are explained in @sec-wave-2D3D-impl, in
+particular @sec-wave2D3D-impl-scalar and
+@sec-wave2D3D-impl-vectorized.
 
 Although vectorization can bring down the CPU time dramatically
 compared with scalar code, there is still some factor 5-10 to win in
@@ -1519,7 +1505,7 @@ to a single index.
 ## The Fortran subroutine
 
 We write a Fortran subroutine `advance` in a file
-[`wave2D_u0_loop_f77.f`](https://github.com/hplgit/fdm-book/tree/master/src/softeng2/wave2D_u0_loop_f77.f)
+[`wave2D_u0_loop_f77.f`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave2D_u0_loop_f77.f)
 for implementing the updating formula
 (@eq-wave-2D3D-impl1-2Du0-ueq-discrete) and setting the solution to zero
 at the boundaries:
@@ -1867,13 +1853,13 @@ void advance(double* u, double* u_1, double* u_2, double* f,
 ## The Cython interface file
 
 All the code above appears in the file
-[`wave2D_u0_loop_c.c`](https://github.com/hplgit/fdm-book/tree/master/src/softeng2/wave2D_u0_loop_c.c).
+[`wave2D_u0_loop_c.c`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave2D_u0_loop_c.c).
 We need to compile this file together with C wrapper code such that
 `advance` can be called from Python. Cython can be used to generate
 appropriate wrapper code.
 The relevant Cython code for interfacing C is
 placed in a file with extension `.pyx`. This file, called
-[`wave2D_u0_loop_c_cy.pyx`](https://github.com/hplgit/fdm-book/tree/master/src/softeng2/wave2D_u0_loop_c_cy.pyx), looks like
+[`wave2D_u0_loop_c_cy.pyx`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave2D_u0_loop_c_cy.pyx), looks like
 
 ```python
 import numpy as np
@@ -1903,7 +1889,7 @@ def advance_cwrap(
 
 We first declare the C functions to be interfaced.
 These must also appear in a C header file,
-[`wave2D_u0_loop_c.h`](https://github.com/hplgit/fdm-book/tree/master/src/softeng2/wave2D_u0_loop_c.h),
+[`wave2D_u0_loop_c.h`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave2D_u0_loop_c.h),
 
 ```python
 extern void advance(double* u, double* u_n, double* u_nm1, double* f,
@@ -2047,8 +2033,8 @@ interfaces for a wide range of
 languages, including Python, Perl, Ruby, and Java.
 However, SWIG is a comprehensive tool with a correspondingly
 steep learning curve. Alternative tools, such as
-[Boost Python](http://www.boost.org/doc/libs/1_51_0/libs/python/doc/index.html), [SIP](http://riverbankcomputing.co.uk/software/sip/intro),
-and [Shiboken](http://qt-project.org/wiki/Category:LanguageBindings::PySide::Shiboken)
+[Boost Python](https://www.boost.org/doc/libs/release/libs/python/doc/html/index.html), [SIP](https://riverbankcomputing.com/software/sip/),
+and [Shiboken](https://wiki.qt.io/Qt_for_Python)
 are similarly comprehensive. Simpler tools include
 [PyBindGen](http://code.google.com/p/pybindgen/).
 

chapters/appendices/softeng2/softeng2.qmd

@@ -401,7 +401,7 @@ D1u + D2u*dt/2 + D3u*dt**2/6 + D4u*dt**3/24
 ```
 The truncation error consists of the terms after the first one ($u'$).
 
-The module file [`trunc/truncation_errors.py`](https://github.com/hplgit/fdm-book/tree/master/src/trunc/truncation_errors.py) contains another class `DiffOp` with symbolic expressions for
+The module file [`trunc/truncation_errors.py`](https://github.com/devitocodes/devito_book/tree/main/src/trunc/truncation_errors.py) contains another class `DiffOp` with symbolic expressions for
 most of the truncation errors listed in the previous section.
 For example:
 
@@ -777,7 +777,7 @@ we pick the rate $r_{m-2}$, which we believe is the best estimation since
 it is based on the two finest meshes.
 
 The `estimate` function is available in a module
-[`trunc_empir.py`](https://github.com/hplgit/fdm-book/tree/master/src/trunc/trunc_empir.py).
+[`trunc_empir.py`](https://github.com/devitocodes/devito_book/tree/main/src/trunc/trunc_empir.py).
 Let us apply this function to estimate the truncation
 error of the Forward Euler scheme. We need a function `decay_FE(dt, N)`
 that can compute (@eq-trunc-decay-FE-R-comp) at the
@@ -814,7 +814,7 @@ The agreement between the theoretical formula (@eq-trunc-decay-FE-R)
 and the computed quantity (ref(@eq-trunc-decay-FE-R-comp)) is
 very good, as illustrated in
 Figures @fig-trunc-fig-FE-rates and @fig-trunc-fig-FE-error.
-The program [`trunc_decay_FE.py`](https://github.com/hplgit/fdm-book/tree/master/src/trunc/trunc_decay_FE.py)
+The program [`trunc_decay_FE.py`](https://github.com/devitocodes/devito_book/tree/main/src/trunc/trunc_decay_FE.py)
 was used to perform the simulations and it can easily be modified to
 test other schemes (see also Exercise @sec-trunc-exer-decay-estimate).
 
@@ -1532,6 +1532,7 @@ in both equations:
 $$
 R_u^{n-\half}= \Oof{\Delta t^2}, \quad R_v^n = \Oof{\Delta t^2}\tp
 $$
+
 ## Linear wave equation in 1D {#sec-trunc-wave-1D}
 
 The standard, linear wave equation in 1D for a function $u(x,t)$ reads
@@ -2234,7 +2235,7 @@ Section @sec-trunc-vib-undamped describes two ways of discretizing
 the initial condition $u'(0)=V$ for a vibration model
 $u''+\omega^2u=0$: a centered difference $[D_{2t}u=V]^0$ or
 a forward difference $[D_t^+u=V]^0$.
-The program [`vib_undamped.py`](https://github.com/hplgit/fdm-book/tree/master/src/vib/vib_undamped.py)
+The program [`vib_undamped.py`](https://github.com/devitocodes/devito_book/tree/main/src/vib/vib_undamped.py)
 solves $u''+\omega^2u=0$ with $[D_{2t}u=0]^0$ and features
 a function `convergence_rates` for computing the order of the
 error in the numerical solution. Modify this program such

chapters/appendices/trunc/trunc.qmd

@@ -49,29 +49,41 @@ This approach has several limitations:
 With Devito, the same problem becomes:
 
 ```python
-from devito import Grid, TimeFunction, Eq, Operator
+from devito import Grid, TimeFunction, Eq, Operator, solve, Constant
+
+# Problem parameters
+Nx = 100
+L = 1.0
+alpha = 1.0  # diffusion coefficient
+F = 0.5      # Fourier number (for stability, F <= 0.5)
+
+# Compute dt from stability condition: F = alpha * dt / dx^2
+dx = L / Nx
+dt = F * dx**2 / alpha
 
 # Create computational grid
-grid = Grid(shape=(101,), extent=(1.0,))
+grid = Grid(shape=(Nx + 1,), extent=(L,))
 
 # Define the unknown field
 u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
 
 # Set initial condition
-u.data[0, 50] = 1.0
+u.data[0, Nx // 2] = 1.0
 
-# Define the PDE update equation
-eq = Eq(u.forward, u + alpha * dt * u.dx2)
+# Define the PDE symbolically and solve for u.forward
+a = Constant(name='a')
+pde = u.dt - a * u.dx2
+update = Eq(u.forward, solve(pde, u.forward))
 
 # Create and run the operator
-op = Operator([eq])
-op(time=1000, dt=dt)
+op = Operator([update])
+op(time=1000, dt=dt, a=alpha)
 ```
 
 This approach offers significant advantages:
 
-1. **Mathematical clarity**: The equation `u.forward = u + alpha * dt * u.dx2`
-   directly mirrors the mathematical formulation
+1. **Mathematical clarity**: The PDE `u.dt - a * u.dx2 = 0` is written symbolically,
+   and Devito derives the update formula automatically using `solve()`
 2. **Automatic optimization**: Devito generates C code with loop tiling,
    SIMD vectorization, and OpenMP parallelization
 3. **Dimension-agnostic**: The same code structure works for 1D, 2D, or 3D

chapters/devito_intro/what_is_devito.qmd

@@ -216,9 +216,7 @@ make Taylor expansions of $A/\Aex$ to see the error more analytically.
 
 As an alternative to examining the accuracy of the damping of a wave
 component, we can perform a general truncation error analysis as
-explained in ref[Appendix @sec-ch-trunc][ in
-[@Langtangen_deqbook_trunc]]["Truncation error analysis": ""
-[@Langtangen_deqbook_trunc]]. Such results are more general, but
+explained in @sec-ch-trunc. Such results are more general, but
 less detailed than what we get from the wave component analysis.  The
 truncation error can almost always be computed and represents the
 error in the numerical model when the exact solution is substituted
@@ -324,17 +322,11 @@ $$
 $$
 ### Truncation error
 We follow the theory explained in
-ref[Appendix @sec-ch-trunc][ in
-[@Langtangen_deqbook_trunc]]["Truncation error analysis": ""
-[@Langtangen_deqbook_trunc]]. The recipe is to set up the
+@sec-ch-trunc. The recipe is to set up the
 scheme in operator notation and use formulas from
-ref[Section @sec-trunc-table][ in
-[@Langtangen_deqbook_trunc]]["Overview of leading-order error terms in finite difference formulas": ""
-[@Langtangen_deqbook_trunc]] to derive an expression for
+@sec-trunc-table to derive an expression for
 the residual. The details are documented in
-ref[Section @sec-trunc-diffu-1D][ in
-[@Langtangen_deqbook_trunc]]["Linear diffusion equation in 1D": ""
-[@Langtangen_deqbook_trunc]]. We end up with a truncation error
+@sec-trunc-diffu-1D. We end up with a truncation error
 $$
 R^n_i = \Oof{\Delta t} + \Oof{\Delta x^2}\tp
 $$
@@ -423,9 +415,7 @@ $F\leq\half$.
 
 ### Truncation error
 The truncation error is derived in
-ref[Section @sec-trunc-diffu-1D][ in
-[@Langtangen_deqbook_trunc]]["Linear diffusion equation in 1D": ""
-[@Langtangen_deqbook_trunc]]:
+@sec-trunc-diffu-1D:
 $$
 R^{n+\half}_i = \Oof{\Delta x^2} + \Oof{\Delta t^2}\tp
 $$

chapters/diffu/diffu_analysis.qmd

@@ -195,7 +195,7 @@ Section @sec-diffu-pde1-analysis.
 
 ## Implementation {#sec-diffu-pde1-FE-code}
 
-The file [`diffu1D_u0.py`](https://github.com/hplgit/fdm-book/tree/master/src/diffu/diffu1D_u0.py)
+The file [`diffu1D_u0.py`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/diffu1D_u0.py)
 contains a complete function `solver_FE_simple`
 for solving the 1D diffusion equation with $u=0$ on the boundary
 as specified in the algorithm above:
@@ -1015,7 +1015,7 @@ and a smooth Gaussian function,
 $$
 I(x) = e^{-\frac{1}{2\sigma^2}(x-L/2)^2}\tp
 $$
-The functions `plug` and `gaussian` in [`diffu1D_u0.py`](https://github.com/hplgit/fdm-book/tree/master/src/diffu/diffu1D_u0.py) run the two cases,
+The functions `plug` and `gaussian` in [`diffu1D_u0.py`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/diffu1D_u0.py) run the two cases,
 respectively:
 
 ```python
@@ -1409,7 +1409,7 @@ The `scipy.sparse.linalg.spsolve` function utilizes the sparse storage
 structure of `A` and performs, in this case, a very efficient Gaussian
 elimination solve.
 
-The program [`diffu1D_u0.py`](https://github.com/hplgit/fdm-book/tree/master/src/diffu/diffu1D_u0.py)
+The program [`diffu1D_u0.py`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/diffu1D_u0.py)
 contains a function `solver_BE`, which implements the Backward Euler scheme
 sketched above.
 As mentioned in Section @sec-diffu-pde1-FE,

chapters/diffu/diffu_fd1.qmd

@@ -108,7 +108,7 @@ def solver_theta(I, a, L, Nx, D, T, theta=0.5, u_L=1, u_R=0,
         u_n, u = u, u_n
 ```
 
-The code is found in the file [`diffu1D_vc.py`](https://github.com/hplgit/fdm-book/tree/master/src/diffu/diffu1D_vc.py).
+The code is found in the file [`diffu1D_vc.py`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/diffu1D_vc.py).
 
 ## Stationary solution {#sec-diffu-varcoeff-stationary}
 
@@ -171,7 +171,7 @@ $x$ lies, and then start evaluating a formula like
 (@eq-diffu-fd2-pde-st-sol-pc). In Python, vectorized expressions may
 help to speed up the computations.
 The convenience classes `PiecewiseConstant` and
-`IntegratedPiecewiseConstant` in the [`Heaviside`](https://github.com/hplgit/fdm-book/tree/master/src/diffu/Heaviside.py)
+`IntegratedPiecewiseConstant` in the [`Heaviside`](https://github.com/devitocodes/devito_book/tree/main/src/diffu/Heaviside.py)
 module were made to simplify programming with
 functions like (@eq-diffu-fd2-pde-st-pc-alpha) and expressions like
 (@eq-diffu-fd2-pde-st-sol-pc). These utilities not only represent