Fixing in docs

demo/unitdisc_helmholtz.py

@@ -104,10 +104,10 @@
 gij = T.coors.get_covariant_metric_tensor()
 ui, vi = TT.local_mesh(True, kind='uniform')
 # Evaluate metric on computational mesh
-g = np.array(sp.lambdify(psi, gij)(ui, vi), dtype=object)
-errorg = inner(1, (du[0]-gradue[0])**2*g[0, 0]+ (du[1]-gradue[1])**2*g[1, 1])
-print('Error gradient = %2.6e'%(np.sqrt(errorg)))
-assert np.sqrt(errorg) < 1e-7
+gij[0, 0] = np.array(sp.lambdify(psi, gij[0, 0])(ui, vi), dtype=object)
+errorg = inner(1, (du[0]-gradue[0])**2*gij[0, 0]+ (du[1]-gradue[1])**2*gij[1, 1])
+print('Error gradient = %2.6e'%(np.sqrt(float(errorg))))
+assert np.sqrt(float(errorg)) < 1e-7
 
 if 'pytest' not in os.environ:
     import matplotlib.pyplot as plt
@@ -133,7 +133,11 @@
     plt.figure()
     plt.contourf(xp, yp, up)
     b = T.coors.get_covariant_basis()
-    bij = np.array(sp.lambdify(psi, b)(ui, vi))
+    #bij = np.array(sp.lambdify(psi, b)(ui, vi)) # no longer allowed in sympy
+    bij = np.zeros((2, 2, N, N))
+    for i in (0, 1):
+        for j in (0, 1):
+            bij[i, j] = sp.lambdify(psi, b[i, j])(ui, vi)
     # Compute Cartesian gradient
     gradu = du[0]*bij[0] + du[1]*bij[1]
     plt.quiver(xi, yi, gradu[0], gradu[1], scale=40, pivot='mid', color='white')

docs/Makefile

@@ -10,7 +10,7 @@ BUILDDIR      = build
 DEMO          = Poisson/poisson.rst \
                 KleinGordon/kleingordon.rst \
                 Poisson3D/poisson3d.rst \
-				PolarHelmholtz/polarhelmholtz.rst \
+				PolarHelmholtz/polarhelmholtz.ipynb \
 				SphereHelmholtz/sphericalhelmholtz.ipynb \
 				KuramatoSivashinsky/kuramatosivashinsky.rst \
 				Stokes/stokes.rst \

docs/demos/DrivenCavity/drivencavity.do.txt

@@ -32,11 +32,11 @@ The nonlinear steady Navier Stokes equations are given in strong form as
 !bt
 \begin{align*}
 \nu \nabla^2 \bs{u} - \nabla p &= \nabla \cdot \bs{u} \bs{u} \quad \text{in }  \Omega , \\
-\nabla \cdot \bs{u} &= 0 \quad \text{in } \Omega  \\
-\int_{\Omega} p dx &= 0 \\
-\bs{u}(x, y=1) = (1, 0) \, &\text{ or }\, \bs{u}(x, y=1) = ((1-x)^2(1+x)^2, 0) \\
-\bs{u}(x, y=-1) &= (0, 0) \\
-\bs{u}(x=\pm 1, y) &= (0, 0)
+\nabla \cdot \bs{u} &= 0 \quad \text{in } \Omega,  \\
+\int_{\Omega} p dx &= 0, \\
+\bs{u}(x, y=1) = (1, 0) \, &\text{ or }\, \bs{u}(x, y=1) = ((1-x)^2(1+x)^2, 0), \\
+\bs{u}(x, y=-1) &= (0, 0), \\
+\bs{u}(x=\pm 1, y) &= (0, 0),
 \end{align*}
 !et
 where $\bs{u}, p$ and $\nu$ are, respectively, the
@@ -141,7 +141,7 @@ $D_1^{N_1}(y)=\text{span}\{\mathcal{Y}_l(y)\}_{l=0}^{N_1-1}$ as
 D1Y = FunctionSpace(N[1], family, quad=quad, bc=(0, 1))
 !ec
 
-where `bc=(1, 0)` fixes the values for $y=1$ and $y=-1$, respectively.
+where `bc=(0, 1)` fixes the values for $y=-1$ and $y=1$, respectively.
 For a regularized lid driven cavity the velocity of the top lid is
 $(1-x)^2(1+x)^2$ and not unity. To implement this boundary condition
 instead, we can make use of "sympy": "https://www.sympy.org" and
@@ -199,7 +199,15 @@ V1 = TensorProductSpace(comm, (D0X, D1Y))
 V0 = TensorProductSpace(comm, (D0X, D0Y))
 P = TensorProductSpace(comm, (PX, PY), modify_spaces_inplace=True)
 !ec
-These tensor product spaces are all scalar valued.
+where ``modify_spaces_inplace=True`` makes use of ``PX`` and
+``PY`` directly. These spaces have now been modified to fit in
+the TensorProductSpace ``P``, along two different directions and as such
+the original spaces have been modified. The default behavior in shenfun
+is to make copies of the 1D spaces under the hood, and thus leave ``PX``
+and ``PY`` untouched. In that case the two new and modified spaces would be accessible
+from ``P.bases``.
+
+Note that all these tensor product spaces are scalar valued.
 The velocity is a vector, and a vector requires a mixed vector basis like
 $W_1^{\bs{N}} = V_1^{\bs{N}} \times V_0^{\bs{N}}$. The vector basis is created
 in shenfun as

docs/demos/KleinGordon/kleingordon.do.txt

@@ -64,7 +64,8 @@ t=0)}{\partial t} = u_t^0. label{eq:init}
 
 The spatial coordinates are here denoted as $\boldsymbol{x} = (x, y, z)$, and
 $t$ is time. The parameter $\gamma=\pm 1$ determines whether the equations are focusing
-($+1$) or defocusing ($-1$) (in the movie we have used $\gamma=1$). The domain $\Omega=[-2\pi, 2\pi]^3$ is triply
+($+1$) or defocusing ($-1$) (in the movie we have used $\gamma=1$).
+The domain $\Omega=[-2\pi, 2\pi]^3$ is triply
 periodic and initial conditions will here be set as
 
 !bt
@@ -113,23 +114,25 @@ the spectral Galerkin method. Being a Galerkin method, we need to reshape the
 governing equations into proper variational forms, and this is done by
 multiplying  (ref{eq:df}) and (ref{eq:du}) with the complex conjugate of proper
 test functions and then integrating
-over the domain. To this end we use testfunctions $g\in V(\Omega)$
-with Eq. (ref{eq:df})  and  $v \in V(\Omega)$ with Eq. (ref{eq:du}), where
-$V(\omega)$ is a suitable function space, and obtain
+over the domain. To this end we make use of the triply periodic tensor product
+function space $W^{\boldsymbol{N}}(\Omega)$ (defined in Eq. (ref{eq:kg:Wn}))
+and use testfunctions $g \in W^{\boldsymbol{N}}$
+with Eq. (ref{eq:df})  and  $v \in W^{\boldsymbol{N}}$ with Eq. (ref{eq:du}),
+and obtain
 
 !bt
 \begin{align}
-\frac{\partial}{\partial t} \int_{\Omega} f\, \overline{g}\, w \,dx &= \int_{\Omega}
-\left(\nabla^2 u - \gamma( u\, - u|u|^2) \right) \overline{g} \, w \,dx, label{eq:df_var} \\
+\frac{\partial}{\partial t} \int_{\Omega} f\, \overline{g}\, w \,d\Omega &= \int_{\Omega}
+\left(\nabla^2 u - \gamma( u\, - u|u|^2) \right) \overline{g} \, w \,d\Omega, label{eq:df_var} \\
 \frac{\partial }{\partial t} \int_{\Omega} u\, \overline{v}\, w \, dx &=
-\int_{\Omega} f\, \overline{v} \, w \, dx. label{eq:kg:du_var}
+\int_{\Omega} f\, \overline{v} \, w \, d\Omega. label{eq:kg:du_var}
 \end{align}
 !et
 Note that the overline is used to indicate a complex conjugate, and
 $w$ is a weight function associated with the test functions. The functions
 $f$ and $u$ are now
 to be considered as trial functions, and the integrals over the
-domain are often referred to as inner products. With inner product notation
+domain are referred to as inner products. With inner product notation
 
 !bt
 \[
@@ -152,7 +155,10 @@ still left open. There are numerous different approaches that one could take for
 discretizing in time, and the first two terms on the right hand side of
 (ref{eq:df_var2}) can easily be treated implicitly as well as explicitly. However,
 the approach we will follow in Sec. (ref{sec:rk}) is a fully explicit 4th order "Runge-Kutta":
-"https://en.wikipedia.org/wiki/Runge-Kutta_methods" method.
+"https://en.wikipedia.org/wiki/Runge-Kutta_methods" method. Also note that
+the inner product in the demo will be computed numerically with quadrature
+through fast Fourier transforms, and the integrals are thus not computed exactly
+for all terms.
 
 ===== Discretization =====
 To find a numerical solution we need to discretize the continuous problem
@@ -193,20 +199,21 @@ for simplicity. A 3D tensor product basis function is now defined as
 \begin{equation}
 \Phi_{lmn}(x,y,z) = e^{\imath \underline{l} x} e^{\imath \underline{m} y}
 e^{\imath \underline{n} z} = e^{\imath
-(\underline{l}x + \underline{m}y + \underline{n}z)}
+(\underline{l}x + \underline{m}y + \underline{n}z)},
 \end{equation}
 !et
 where the indices for $y$- and $z$-direction are $\underline{m}=\frac{2\pi}{L}m,
 \underline{n}=\frac{2\pi}{L}n$, and $\boldsymbol{m}$ and $\boldsymbol{n}$ are the same as
 $\boldsymbol{l}$ due to using the same number of basis functions for each direction. One
 distinction, though, is that for the $z$-direction expansion coefficients are only stored for
-$n=0, 1, \ldots, N/2$ due to Hermitian symmetry (real input data).
+$n=0, 1, \ldots, N/2$ due to Hermitian symmetry (real input data). However, for simplicity,
+we still write the sum in Eq. (ref{eq:usg}) over the entire range of basis functions.
 
 We now look for solutions of the form
 
 !bt
 \begin{equation}
-u(x, y, z, t) = \sum_{l=-N/2}^{N/2-1}\sum_{m=-N/2}^{N/2-1}\sum_{n=0}^{N/2}
+u(x, y, z, t) = \sum_{l=-N/2}^{N/2-1}\sum_{m=-N/2}^{N/2-1}\sum_{n=-N/2}^{N/2-1}
 \hat{u}_{lmn} (t)\Phi_{lmn}(x,y,z). label{eq:usg}
 \end{equation}
 !et
@@ -232,13 +239,13 @@ is set to $[-2\pi, 2\pi]^3$ and not the more common $[0, 2\pi]^3$. We have
 
 !bt
 \begin{equation}
-\boldsymbol{u} = \mathcal{F}_z^{-1}\left(\mathcal{F}_y^{-1}\left(\mathcal{F}_z^{-1}\left(\hat{\boldsymbol{u}}\right)\right)\right) label{eq:uxyz}
+\boldsymbol{u} = \mathcal{F}_x^{-1}\left(\mathcal{F}_y^{-1}\left(\mathcal{F}_z^{-1}\left(\hat{\boldsymbol{u}}\right)\right)\right) label{eq:uxyz}
 \end{equation}
 !et
 
 with $\boldsymbol{u} = \{u(x_i, y_j, z_k)\}_{(i,j,k)\in \boldsymbol{i} \times \boldsymbol{j} \times \boldsymbol{k}}$
 and where $\mathcal{F}_x^{-1}$ is the inverse Fourier transform along the direction $x$, for
-all indices in the other direction, i.e., for $(j, k) \in \boldsymbol{j} \times \boldsymbol{k}$. Note that the three
+all indices in the other direction. Note that the three
 inverse FFTs are performed sequentially, one direction at the time, and that there is no
 scaling factor due to
 the definition used for the inverse "Fourier transform": "https://mpi4py-fft.readthedocs.io/en/latest/dft.html"
@@ -270,14 +277,14 @@ unity:
 \int_{0}^{L} e^{\imath \underline{k}x} e^{- \imath \underline{l}x} \frac{1}{L} dx = \delta_{kl}.
 !et
 
-With this weight function the scalar product and the forward transform
+With this weight function the (discrete) scalar product and the forward transform
 are the same and we obtain:
 
 !bt
 \begin{equation}
 \left(u, v \right) = \boldsymbol{\hat{u}} =
 \left(\frac{1}{N}\right)^3
-\mathcal{F}_z\left(\mathcal{F}_y\left(\mathcal{F}_x\left(\boldsymbol{u}\right)\right)\right),
+\mathcal{F}_z\left(\mathcal{F}_y\left(\mathcal{F}_x\left(\boldsymbol{u}\right)\right)\right).
 \end{equation}
 !et
 
@@ -322,8 +329,8 @@ for all $t>0$, find $(f, u) \in W^N \times W^N$  such that
 !bt
 \begin{align}
 \frac{\partial}{\partial t} (f, g) &= -(\nabla u, \nabla g)
--\gamma \left( u - u|u|^2, g \right), label{eq:dff} \\
-\frac{\partial }{\partial t} (u, v) &= (f, v) \quad \forall \, (g, v) \in W^N \times W^N. label{eq:kg:duu}
+-\gamma \left( u - u|u|^2, g \right) \quad \forall \, g \in W^{N}, label{eq:dff} \\
+\frac{\partial }{\partial t} (u, v) &= (f, v) \quad \forall \, v \in W^N. label{eq:kg:duu}
 \end{align}
 !et
 where $u(x, y, z, 0)$ and $f(x, y, z, 0)$ are given as the initial conditions

docs/demos/KuramatoSivashinsky/kuramatosivashinsky.do.txt

@@ -55,8 +55,8 @@ the spectral Galerkin method. Being a Galerkin method, we need to reshape the
 governing equations into proper variational forms, and this is done by
 multiplying  (ref{eq:ks}) with the complex conjugate of a proper
 test function and then integrating
-over the domain. To this end we use testfunction $v\in V(\Omega)$, where $V(\Omega)$
-is some suitable function space, and obtain
+over the domain. To this end we use testfunction $v\in W^N(\Omega)$, where $W^N(\Omega)$
+is a suitable function space (defined in Eq. (ref{eq:Wn})), and obtain
 
 !bt
 \begin{equation}
@@ -159,13 +159,14 @@ is set to $[-30\pi, 30\pi]^2$ and not the more common $[0, 2\pi]^2$. We now have
 
 !bt
 \begin{equation}
-u(x_i, y_j) =
-\mathcal{F}_y^{-1}\left(\mathcal{F}_x^{-1}\left(\hat{u}\right)\right)
-\, \forall\, (i,j)\in\boldsymbol{i} \times \boldsymbol{j},
+\boldsymbol{u} =
+\mathcal{F}_x^{-1}\left(\mathcal{F}_y^{-1}\left(\boldsymbol{\hat{u}}\right)\right),
 \end{equation}
 !et
-where $\mathcal{F}_x^{-1}$ is the inverse Fourier transform along direction
-$x$, for all $j \in \boldsymbol{j} $. Note that the two
+where $\boldsymbol{u} = \{u(x_i, y_j)\}_{(i,j)\in \boldsymbol{i} \times \boldsymbol{j}}$,
+$\boldsymbol{\hat{u}} = \{\hat{u}_{lm}\}_{(l,m)\in \boldsymbol{l} \times \boldsymbol{m}}$
+and $\mathcal{F}_x^{-1}$ is the inverse Fourier transform along direction
+$x$, for all indices in the other direction. Note that the two
 inverse FFTs are performed sequentially, one direction at the time, and that
 there is no scaling factor due
 the definition used for the inverse
@@ -186,10 +187,9 @@ computed using forward FFTs (using weight functions $w=1/L$):
 
 !bt
 \begin{equation}
-\left(u, \Phi_{lm}\right) = \hat{u}_{lm} =
+\boldsymbol{\hat{u}} =
 \frac{1}{N^2}
-\mathcal{F}_l\left(\mathcal{F}_m\left({u}\right)\right)
-\quad \forall (l,m) \in \boldsymbol{l} \times \boldsymbol{m},
+\mathcal{F}_y\left(\mathcal{F}_x\left(\boldsymbol{u}\right)\right),
 \end{equation}
 !et
 
@@ -200,13 +200,14 @@ exact derivatives of the nabla operator, we have
 !bt
 \begin{align}
 (\nabla^2 u, v) &=
--(\underline{l}^2+\underline{m}^2)\hat{u}_{lm}, \\
-(\nabla^4 u, v) &= (\underline{l}^2+\underline{m}^2)^2\hat{u}_{lm}, \\
-(|\nabla u|^2, v) &= \widehat{|\nabla u|^2}_{lm}
+\left(-(\underline{l}^2+\underline{m}^2)\hat{u}_{lm}\right), \\
+(\nabla^4 u, v) &= \left((\underline{l}^2+\underline{m}^2)^2\hat{u}_{lm}\right), \\
+(|\nabla u|^2, v) &= \left(\widehat{|\nabla u|^2}_{lm}\right),
 \end{align}
 !et
 
-and as such the equation to be solved for each wavenumber can be found directly as
+where the indices on the right hand side run over $\boldsymbol{l} \times \boldsymbol{m}$.
+We find that the equation to be solved for each wavenumber can be found directly as
 
 !bt
 \begin{equation}
@@ -281,8 +282,8 @@ with the second when computing odd derivatives.
 
 !bt
 \begin{align}
-u & = \sum_{k=-N/2}^{N/2-1} \hat{u} e^{\imath k x}\\
-u & = \sideset{}{'}\sum_{k=-N/2}^{N/2} \hat{u} e^{\imath k x}
+u & = \sum_{k=-N/2}^{N/2-1} \hat{u}_k e^{\imath k x}\\
+u & = \sideset{}{'}\sum_{k=-N/2}^{N/2} \hat{u}_k e^{\imath k x}
 \end{align}
 !et
 Here $\sideset{}{'}\sum$ means that the first and last items in the sum are

docs/demos/Poisson3D/poisson3d.do.txt

@@ -35,9 +35,8 @@ u(x, y, 2\pi) &= u(x, y, 0),
 where $u(\bs{x})$ is the solution and $f(\bs{x})$ is a function. The domain
 $\Omega = [-1, 1]\times [0, 2\pi]^2$.
 
-To solve Eq. (ref{eq:3d:poisson}) with the Galerkin method we need smooth basis
-functions, $v(\bs{x})$, that live
-in the Hilbert space $H^1(\Omega)$ and that satisfy the given boundary
+To solve Eq. (ref{eq:3d:poisson}) with the Galerkin method we will make use of
+smooth basis functions, $v(\bs{x})$, that satisfy the given boundary
 conditions. To this end we will use one basis function for the $x$-direction,
 $\mathcal{X}(x)$,
 one for the $y$-direction, $\mathcal{Y}(y)$, and one for the $z$-direction,
@@ -97,7 +96,7 @@ u(\bs{x}) &= \sum_{l\in \bs{l}^{N_0}} \sum_{m\in \bs{m}^{N_1}}\sum_{n\in
 !et
 
 where $\hat{u}_{lmn}$ are components of the expansion coefficients for $u$ and
-the second form, $(\hat{u}_{\ts{k}})_{\ts{k}\in\bs{k}}$, is a shorter,
+the second form, $\{\hat{u}_{\ts{k}}\}_{\ts{k}\in\bs{k}}$, is a shorter,
 simplified notation, with sans-serif $\ts{k}=(l, m, n)$.
 The expansion coefficients are the unknowns in the spectral Galerkin method.
 
@@ -119,8 +118,10 @@ The weighted integrals, weighted by $w(\bs{x})$, are called inner products, and
 \end{equation}
 !et
 The integral can either be computed exactly, or with quadrature. The advantage
-of the latter is that it is generally faster, and that non-linear terms may be
-computed just as quickly as linear. For a linear problem, it does not make much of a difference, if any at all. Approximating the integral with quadrature, we obtain
+of the latter is that it is faster (through Fast Fourier transforms),
+and that non-linear terms may be computed just as quickly as linear.
+For a linear problem, it does not make much of a difference, if any at all.
+Approximating the integral with quadrature, we obtain
 !bt
 \begin{align*}
 \int_{\Omega} u \, \overline{v} \, w\, \bs{dx} &\approx \langle u, v
@@ -129,11 +130,10 @@ computed just as quickly as linear. For a linear problem, it does not make much
 \end{align*}
 !et
 where $w(\bs{x})$ now are the quadrature weights. The quadrature points
-$(x_i)_{i=0}^{N_0-1}$ are specific to the chosen basis, and even within basis there
+$\{x_i\}_{i=0}^{N_0-1}$ are specific to the chosen basis, and even within basis there
 are two different choices based on which quadrature rule is selected, either
-Gauss or Gauss-Lobatto. The quadrature points for the Fourier bases are the
-uniform $(y_j)_{j=0}^{N_1-1}=2\pi j / N_1$ and $(z_k)_{k=0}^{N_2-1} = 2 \pi
-k/N_2$.
+Gauss or Gauss-Lobatto. The quadrature points for the Fourier bases are simply
+uniformly distributed throughout the domain.
 
 Inserting for test function (ref{eq:3d:u}) and trialfunction
 $v_{pqr} = \mathcal{X}_{p} \mathcal{Y}_q \mathcal{Z}_r$ on the
@@ -157,7 +157,7 @@ b_{pl} = \left( \mathcal{X}_l, \mathcal{X}_p \right)_w^{N_0} = \sum_{i=0}^{N_0-1
 \mathcal{X}_p(x_i) w(x_i),
 \end{equation}
 !et
-is used to represent an $L_2$ inner product along only the first, nonperiodic,
+is used to represent a discrete $L_2$ inner product along only the first, nonperiodic,
 direction. The delta functions above come from integrating over the two periodic
 directions, where we use constant weight functions $w=1/(2\pi)$ in the
 inner products
@@ -316,8 +316,8 @@ direction to be local to the processor. If the `SD` basis is the last to be
 transformed, then the data will be aligned in this direction, whereas the other
 two directions may both, or just one of them, be distributed.
 
-Note that `X` is a list containing local values of the arrays $(x_i)_{i=0}^{N_0-1}$,
-$(y_j)_{j=0}^{N_1-1}$ and $(z_k)_{k=0}^{N_2-1}$.
+Note that `X` is a list containing local values of the arrays $\{x_i\}_{i=0}^{N_0-1}$,
+$\{y_j\}_{j=0}^{N_1-1}$ and $\{z_k\}_{k=0}^{N_2-1}$.
 Now, it's not possible to run a jupyter notebook with more than one process,
 but we can imagine running "the complete solver": "https://github.com/spectralDNS/shenfun/blob/master/demo/poisson3D.py"
 with 4 procesors and a processor mesh of shape $2\times 2$.

docs/demos/PolarHelmholtz/polarhelmholtz.do.txt

@@ -78,7 +78,7 @@ A discrete Fourier approximation space with $N$ basis functions is then
 
 !bt
 \begin{equation}
-V_F^N = \text{span} \{\exp(\imath k \theta)\}, \text{ for } k \in K,
+V_F^N = \text{span} \{\exp(\imath k \theta)| \text{ for } k \in K\},
 \end{equation}
 !et
 
@@ -87,7 +87,7 @@ is real, there is Hermitian symmetry and $\hat{u}_{k,j} = \hat{u}_{k,-j}^*$
 (with $*$ denoting a complex conjugate).
 For this reason we use only $k \in K=\{0, 1, \ldots, N/2\}$ in solving for
 $\hat{u}_{kj}$, and then use Hermitian symmetry to get the remaining
-unknowns.
+unknowns. This is handled under the hood by fast Fourier transforms.
 
 The radial basis is more tricky, because there is a nontrivial 'boundary'
 condition (pole condition) that needs to be applied at the center of the disc $(r=0)$
@@ -427,7 +427,10 @@ on the computational grid
 !bc pycod
 ui, vi = TT.local_mesh(True)
 b = T.coors.get_covariant_basis()
-bij = np.array(sp.lambdify(psi, b)(ui, vi))
+bij = np.zeros((2, 2, N, N))
+for i in (0, 1):
+    for j in (0, 1):
+        bij[i, j] = sp.lambdify(psi, b[i, j])(ui, vi)
 gradu = du[0]*bij[0] + du[1]*bij[1]
 !ec
 
@@ -449,10 +452,10 @@ gradue = Array(V, buffer=grad(u).tosympy(basis=ue, psi=psi))
 gij = T.coors.get_covariant_metric_tensor()
 ui, vi = TT.local_mesh(True, kind='uniform')
 # Evaluate metric on computational mesh
-g = np.array(sp.lambdify(psi, gij)(ui, vi), dtype=object)
+gij[0, 0] = sp.lambdify(psi, gij[0, 0])(ui, vi)
 # Compute L2 error
-errorg = inner(1, (du[0]-gradue[0])**2*g[0, 0]+ (du[1]-gradue[1])**2*g[1, 1])
-print('Error gradient', np.sqrt(errorg))
+errorg = inner(1, (du[0]-gradue[0])**2*gij[0, 0]+ (du[1]-gradue[1])**2*gij[1, 1])
+print('Error gradient', np.sqrt(float(errorg)))
 !ec
 
 % if FORMAT not in ("ipynb"):

docs/demos/RayleighBenard/rayleighbenard.do.txt

@@ -65,7 +65,7 @@ $\bs{H} = (H_x, H_y) = (\bs{u} \cdot \nabla) \bs{u}$
 \end{equation}
 !et
 
-This equation is solved with $u(\pm 1) = \partial u/\partial x(\pm 1) = 0$, where the latter follows from the
+This equation is solved with $u(\pm 1,y,t) = \partial u/\partial x(\pm 1,y,t) = 0$, where the latter follows from the
 divergence constraint. In summary, we now seem to have the following equations to solve:
 
 !bt