Clean up example taus, modify upgrade docs

docs/pages/installation.rst

@@ -79,20 +79,20 @@ Alternatively, you can clone the d3 branch from the source repository and instal
 
     git clone -b d3 https://github.com/DedalusProject/dedalus
     cd dedalus
-    pip3 install .
+    pip3 install --no-cache .
 
 Updating Dedalus
 ----------------
 
 If Dedalus was installed using the conda script or from GitHub with pip, it can also be updated using pip::
 
-    pip3 install --upgrade --no-cache http://github.com/dedalusproject/dedalus/zipball/d3/
+    pip3 install --upgrade --force-reinstall --no-deps --no-cache http://github.com/dedalusproject/dedalus/zipball/d3/
 
 If Dedalus was built from a clone of the source repository, first pull new changes and then reinstall with pip::
 
     cd /path/to/dedalus/repo
     git pull
-    pip3 install --upgrade --force-reinstall .
+    pip3 install --upgrade --force-reinstall --no-deps --no-cache .
 
 **Note**: any custom FFTW/MPI paths set in the conda script or during the original installation will also need to be exported for the upgrade commands to work.
 

docs/pages/tau_method.rst

@@ -121,9 +121,9 @@ We will also include a constant scalar tau term that will allow us to impose the
     # Fields
     p = dist.Field(name='p', bases=(xbasis,ybasis))
     u = dist.VectorField(coords, name='u', bases=(xbasis,ybasis))
-    tau_1 = dist.VectorField(coords, name='tau_1', bases=xbasis)
-    tau_2 = dist.VectorField(coords, name='tau_2', bases=xbasis)
-    tau_c = dist.Field(name='tau_c')
+    tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
+    tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)
+    tau_p = dist.Field(name='tau_p')
 
 We then create substitutions for :math:`G` and :math:`P(y)`.
 Specification of and multiplication by :math:`P(y)` are handled through the ``Lift`` operator, which here simply multiplies its argument by the specified mode/element of a selected basis.
@@ -135,17 +135,17 @@ Here we'll take :math:`P(y)` to be the highest mode in the Chebyshev-U basis, in
     ex, ey = coords.unit_vector_fields(dist)
     lift_basis = ybasis.clone_with(a=1/2, b=1/2) # Chebyshev U basis
     lift = lambda A, n: d3.Lift(A, lift_basis, -1) # Shortcut for multiplying by U_{N-1}(y)
-    grad_u = d3.grad(u) - ey*lift(tau_1) # Operator representing G
+    grad_u = d3.grad(u) - ey*lift(tau_u1) # Operator representing G
 
 We can then create a problem and enter the PDE, boundary condtions, and pressure gauge in vectorial form using these substitutions.
 Note here we will add the contant tau term to the divergence equation, which introduces a degree of freedom allowing the imposition of the pressure gauge (otherwise the integral of the divergence equation will be redundant with integrals of the inflow boundary conditions).
 
 .. code-block:: python
 
     # Problem
-    problem = d3.IVP([p, u, tau_1, tau_2, tau_c], namespace=locals())
-    problem.add_equation("trace(grad_u) + tau_c = 0")
-    problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) + lift(tau_2) = f")
+    problem = d3.IVP([p, u, tau_u1, tau_u2, tau_p], namespace=locals())
+    problem.add_equation("trace(grad_u) + tau_p = 0")
+    problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) + lift(tau_u2) = f")
     problem.add_equation("u(y=-1) = 0")
     problem.add_equation("u(y=+1) = 0")
     problem.add_equation("integ(p) = 0")
@@ -166,28 +166,27 @@ For instance, to enter the above equation set with homogeneous Dirichlet boundar
     # Fields
     p = dist.Field(name='p', bases=disk_basis)
     u = dist.VectorField(coords, name='u', bases=disk_basis)
-    tau_1 = dist.VectorField(coords, name='tau_1', bases=phi_basis)
-    tau_c = dist.Field(name='tau_c')
+    tau_u = dist.VectorField(coords, name='tau_u', bases=phi_basis)
+    tau_p = dist.Field(name='tau_p')
 
 The disk and ball bases are not direct-product bases, so the tau terms can't actually be written just as the tau field times a radial polynomial.
 Instead, for each horizontal mode (azimuthal mode :math:`m` in the disk and spherical harmonic :math:`\ell` in the ball), the tau term is multiplied by the highest degree radial polynomial in the basis for that particular mode.
 The ``Lift`` operator does this under the hood, and is why we use it rather than explicitly writing out the tau polynomials.
-Here, we'll use a tau polynomial from the second-derivative basis:
+We've found that using tau polynomials from the original bases seems to give good results in the disk and ball:
 
 .. code-block:: python
 
     # Substitutions
-    lift_basis = disk_basis.clone_with(k=2) # Second-derivative basis
-    lift = lambda A, n: d3.Lift(A, lift_basis, -1)
+    lift = lambda A, n: d3.Lift(A, disk_basis, -1)
 
 Now we can enter the PDE with just the single tau term in the momentum equation:
 
 .. code-block:: python
 
     # Problem
-    problem = d3.IVP([p, u, tau_1, tau_c], namespace=locals())
-    problem.add_equation("div(u) + tau_c = 0")
-    problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_1) = f")
+    problem = d3.IVP([p, u, tau_u, tau_p], namespace=locals())
+    problem.add_equation("div(u) + tau_p = 0")
+    problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = f")
     problem.add_equation("u(r=1) = 0")
     problem.add_equation("integ(p) = 0")
 

examples/evp_1d_waves_on_a_string/waves_on_a_string.py

@@ -35,19 +35,19 @@
 
 # Fields
 u = dist.Field(name='u', bases=xbasis)
-tau1 = dist.Field(name='tau1')
-tau2 = dist.Field(name='tau2')
+tau_1 = dist.Field(name='tau_1')
+tau_2 = dist.Field(name='tau_2')
 s = dist.Field(name='s')
 
 # Substitutions
 dx = lambda A: d3.Differentiate(A, xcoord)
 lift_basis = xbasis.clone_with(a=1/2, b=1/2) # First derivative basis
-lift = lambda A, n: d3.Lift(A, lift_basis, n)
-ux = dx(u) + lift(tau1,-1) # First-order reduction
+lift = lambda A: d3.Lift(A, lift_basis, -1)
+ux = dx(u) + lift(tau_1) # First-order reduction
 
 # Problem
-problem = d3.EVP([u, tau1, tau2], eigenvalue=s, namespace=locals())
-problem.add_equation("s*u + dx(ux) + lift(tau2,-1) = 0")
+problem = d3.EVP([u, tau_1, tau_2], eigenvalue=s, namespace=locals())
+problem.add_equation("s*u + dx(ux) + lift(tau_2) = 0")
 problem.add_equation("u(x=0) = 0")
 problem.add_equation("u(x=Lx) = 0")
 

examples/ivp_2d_rayleigh_benard/rayleigh_benard.py

@@ -52,11 +52,11 @@
 p = dist.Field(name='p', bases=(xbasis,zbasis))
 b = dist.Field(name='b', bases=(xbasis,zbasis))
 u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))
-taup = dist.Field(name='taup')
-tau1b = dist.Field(name='tau1b', bases=xbasis)
-tau2b = dist.Field(name='tau2b', bases=xbasis)
-tau1u = dist.VectorField(coords, name='tau1u', bases=xbasis)
-tau2u = dist.VectorField(coords, name='tau2u', bases=xbasis)
+tau_p = dist.Field(name='tau_p')
+tau_b1 = dist.Field(name='tau_b1', bases=xbasis)
+tau_b2 = dist.Field(name='tau_b2', bases=xbasis)
+tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis)
+tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis)
 
 # Substitutions
 kappa = (Rayleigh * Prandtl)**(-1/2)
@@ -65,17 +65,17 @@
 ex, ez = coords.unit_vector_fields(dist)
 integ = lambda A: d3.Integrate(d3.Integrate(A, 'x'), 'z')
 lift_basis = zbasis.clone_with(a=1/2, b=1/2) # First derivative basis
-lift = lambda A, n: d3.Lift(A, lift_basis, n)
-grad_u = d3.grad(u) + ez*lift(tau1u,-1) # First-order reduction
-grad_b = d3.grad(b) + ez*lift(tau1b,-1) # First-order reduction
+lift = lambda A: d3.Lift(A, lift_basis, -1)
+grad_u = d3.grad(u) + ez*lift(tau_u1) # First-order reduction
+grad_b = d3.grad(b) + ez*lift(tau_b1) # First-order reduction
 
 # Problem
 # First-order form: "div(f)" becomes "trace(grad_f)"
 # First-order form: "lap(f)" becomes "div(grad_f)"
-problem = d3.IVP([p, b, u, taup, tau1b, tau2b, tau1u, tau2u], namespace=locals())
-problem.add_equation("trace(grad_u) + taup = 0")
-problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau2b,-1) = - dot(u,grad(b))")
-problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) + lift(tau2u,-1) - b*ez = - dot(u,grad(u))")
+problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
+problem.add_equation("trace(grad_u) + tau_p = 0")
+problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - dot(u,grad(b))")
+problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - dot(u,grad(u))")
 problem.add_equation("b(z=0) = Lz")
 problem.add_equation("u(z=0) = 0")
 problem.add_equation("b(z=Lz) = 0")

examples/ivp_ball_internally_heated_convection/internally_heated_convection.py

@@ -2,9 +2,8 @@
 Dedalus script simulating internally-heated Boussinesq convection in the ball.
 This script demonstrates soving an initial value problem in the ball. It can be
 ran serially or in parallel, and uses the built-in analysis framework to save
-data snapshots to HDF5 files. The `plot_equator.py` and `plot_meridian.py` scripts
-can be used to produce plots from the saved data. The simulation should take
-roughly 15 cpu-minutes to run.
+data snapshots to HDF5 files. The `plot_ball.py` script can be used to produce
+plots from the saved data. The simulation should take roughly 15 cpu-minutes to run.
 
 The strength of gravity is proportional to radius, as for a constant density ball.
 The problem is non-dimensionalized using the ball radius and freefall time, so
@@ -18,8 +17,8 @@
 boundary. The convection is driven by the internal heating term with a conductive
 equilibrium of T(r) = 1 - r**2.
 
-For incompressible hydro in the ball, we need one tau terms for each the velocity
-and temperature. Here we choose to lift them to the natural output (k=2) basis.
+For incompressible hydro in the ball, we need one tau term each for the velocity
+and temperature. Here we choose to lift them to the original (k=0) basis.
 
 The simulation will run to t=10, about the time for the first convective plumes
 to hit the top boundary. After running this initial simulation, you can restart
@@ -28,8 +27,7 @@
 To run, restart, and plot using e.g. 4 processes:
     $ mpiexec -n 4 python3 internally_heated_convection.py
     $ mpiexec -n 4 python3 internally_heated_convection.py --restart
-    $ mpiexec -n 4 python3 plot_equator.py slices/*.h5
-    $ mpiexec -n 4 python3 plot_meridian.py slices/*.h5
+    $ mpiexec -n 4 python3 plot_ball.py slices/*.h5
 """
 
 import sys
@@ -74,8 +72,7 @@
 T_source = 6
 kappa = (Rayleigh * Prandtl)**(-1/2)
 nu = (Rayleigh / Prandtl)**(-1/2)
-lift_basis = basis.clone_with(k=2) # Natural output
-lift = lambda A, n: d3.Lift(A, lift_basis, n)
+lift = lambda A: d3.Lift(A, basis, -1)
 strain_rate = d3.grad(u) + d3.trans(d3.grad(u))
 shear_stress = d3.angular(d3.radial(strain_rate(r=1), index=1))
 integ = lambda A: d3.Integrate(A, coords)
@@ -84,8 +81,8 @@
 # Problem
 problem = d3.IVP([p, u, T, tau_p, tau_u, tau_T], namespace=locals())
 problem.add_equation("div(u) + tau_p = 0")
-problem.add_equation("dt(u) - nu*lap(u) + grad(p) - r_vec*T + lift(tau_u,-1) = - cross(curl(u),u)")
-problem.add_equation("dt(T) - kappa*lap(T) + lift(tau_T,-1) = - dot(u,grad(T)) + kappa*T_source")
+problem.add_equation("dt(u) - nu*lap(u) + grad(p) - r_vec*T + lift(tau_u) = - cross(curl(u),u)")
+problem.add_equation("dt(T) - kappa*lap(T) + lift(tau_T) = - dot(u,grad(T)) + kappa*T_source")
 problem.add_equation("shear_stress = 0")  # Stress free
 problem.add_equation("radial(u(r=1)) = 0")  # No penetration
 problem.add_equation("rad(grad(T)(r=1)) = -2")

examples/ivp_ball_internally_heated_convection/plot_ball.py

@@ -5,7 +5,7 @@
     plot_sphere.py <files>... [--output=<dir>]
 
 Options:
-    --output=<dir>  Output directory [default: ./frames_sphere]
+    --output=<dir>  Output directory [default: ./frames]
 
 """
 

examples/ivp_disk_libration/libration.py

@@ -13,7 +13,7 @@
     nu = Ekman
 
 For incompressible hydro in the disk, we need one tau term for the velocity.
-Here we lift to the natural output (k=2) basis.
+Here we lift to the original (k=0) basis.
 
 To run and plot using e.g. 4 processes:
     $ mpiexec -n 4 python3 libration.py
@@ -54,8 +54,7 @@
 phi, r = dist.local_grids(basis)
 nu = Ekman
 integ = lambda A: d3.Integrate(A, coords)
-lift_basis = basis.clone_with(k=2) # Natural output basis
-lift = lambda A, n: d3.Lift(A, lift_basis, n)
+lift = lambda A: d3.Lift(A, basis, -1)
 
 # Background librating flow
 u0_real = dist.VectorField(coords, bases=basis)
@@ -68,7 +67,7 @@
 # Problem
 problem = d3.IVP([p, u, tau_u, tau_p], time=t, namespace=locals())
 problem.add_equation("div(u) + tau_p = 0")
-problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u,-1) = - dot(u, grad(u0)) - dot(u0, grad(u))")
+problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - dot(u, grad(u0)) - dot(u0, grad(u))")
 problem.add_equation("u(r=1) = 0")
 problem.add_equation("integ(p) = 0")
 

examples/ivp_shell_convection/plot_shell.py

@@ -1,11 +1,11 @@
 """
-Plot spherical snapshots.
+Plot cutaway spherical shell outputs.
 
 Usage:
-    plot_sphere.py <files>... [--output=<dir>]
+    plot_shell.py <files>... [--output=<dir>]
 
 Options:
-    --output=<dir>  Output directory [default: ./frames_sphere]
+    --output=<dir>  Output directory [default: ./frames]
 
 """
 

examples/ivp_shell_convection/shell_convection.py

@@ -2,7 +2,7 @@
 Dedalus script simulating Boussinesq convection in a spherical shell. This script
 demonstrates solving an initial value problem in the shell. It can be ran serially
 or in parallel, and uses the built-in analysis framework to save data snapshots
-to HDF5 files. The `plot_sphere.py` script can be used to produce plots from the
+to HDF5 files. The `plot_shell.py` script can be used to produce plots from the
 saved data. The simulation should take about 10 cpu-minutes to run.
 
 The problem is non-dimensionalized using the shell thickness and freefall time, so
@@ -66,16 +66,16 @@
 rvec = dist.VectorField(coords, bases=basis.radial_basis)
 rvec['g'][2] = r
 lift_basis = basis.clone_with(k=1) # First derivative basis
-lift = lambda A, n: d3.Lift(A, lift_basis, n)
-grad_u = d3.grad(u) + rvec*lift(tau_u1,-1) # First-order reduction
-grad_b = d3.grad(b) + rvec*lift(tau_b1,-1) # First-order reduction
+lift = lambda A: d3.Lift(A, lift_basis, -1)
+grad_u = d3.grad(u) + rvec*lift(tau_u1) # First-order reduction
+grad_b = d3.grad(b) + rvec*lift(tau_b1) # First-order reduction
 integ = lambda A: d3.Integrate(A, coords)
 
 # Problem
 problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals())
 problem.add_equation("trace(grad_u) + tau_p = 0")
-problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2,-1) = - dot(u,grad(b))")
-problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*er + lift(tau_u2,-1) = - dot(u,grad(u))")
+problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - dot(u,grad(b))")
+problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*er + lift(tau_u2) = - dot(u,grad(u))")
 problem.add_equation("b(r=Ri) = 1")
 problem.add_equation("u(r=Ri) = 0")
 problem.add_equation("b(r=Ro) = 0")

examples/ivp_sphere_shallow_water/plot_sphere.py

@@ -1,11 +1,11 @@
 """
-Plot planes from joint analysis files.
+Plot cutaway ball outputs.
 
 Usage:
-    plot_sphere.py <files>... [--output=<dir>]
+    plot_ball.py <files>... [--output=<dir>]
 
 Options:
-    --output=<dir>  Output directory [default: ./frames_sphere]
+    --output=<dir>  Output directory [default: ./frames]
 
 """
 

examples/lbvp_2d_poisson/poisson.py

@@ -32,8 +32,8 @@
 
 # Fields
 u = dist.Field(name='u', bases=(xbasis, ybasis))
-tau1 = dist.Field(name='tau1', bases=xbasis)
-tau2 = dist.Field(name='tau2', bases=xbasis)
+tau_1 = dist.Field(name='tau_1', bases=xbasis)
+tau_2 = dist.Field(name='tau_2', bases=xbasis)
 
 # Forcing
 x, y = dist.local_grids(xbasis, ybasis)
@@ -51,8 +51,8 @@
 lift = lambda A, n: d3.Lift(A, lift_basis, n)
 
 # Problem
-problem = d3.LBVP([u, tau1, tau2], namespace=locals())
-problem.add_equation("lap(u) + lift(tau1,-1) + lift(tau2,-2) = f")
+problem = d3.LBVP([u, tau_1, tau_2], namespace=locals())
+problem.add_equation("lap(u) + lift(tau_1,-1) + lift(tau_2,-2) = f")
 problem.add_equation("u(y=0) = g")
 problem.add_equation("dy(u)(y=Ly) = h")
 

examples/nlbvp_ball_lane_emden/lane_emden.py

@@ -28,7 +28,7 @@
 and R can then be recovered from f(r=0) = R**(2/(n-1)).
 
 For a scalar Laplacian in the ball, we need a single tau term. Here we choose
-to lift it to the natural output (k=2) basis.
+to lift it to the original (k=0) basis.
 
 References:
     [1]: http://en.wikipedia.org/wiki/Lane–Emden_equation
@@ -63,12 +63,11 @@
 tau = dist.Field(name='tau', bases=basis.S2_basis(radius=1))
 
 # Substitutions
-lift_basis = basis.clone_with(k=2) # Natural output basis
-lift = lambda A, n: d3.Lift(A, lift_basis, n)
+lift = lambda A: d3.Lift(A, basis, -1)
 
 # Problem
 problem = d3.NLBVP([f, tau], namespace=locals())
-problem.add_equation("lap(f) + lift(tau,-1) = - f**n")
+problem.add_equation("lap(f) + lift(tau) = - f**n")
 problem.add_equation("f(r=1) = 0")
 
 # Initial guess