Add gauge constraint docs

README.md

@@ -19,20 +19,20 @@
 **Note: this branch is the development head of v3 of the codebase, which is currently under beta-release.
 The development head of v2 of the codebase is on the ["v2_master" branch](https://github.com/DedalusProject/dedalus/tree/v2_master).**
 
-Dedalus is a flexible framework for solving partial differential equations using spectral methods.
+Dedalus is a flexible framework for solving partial differential equations using modern spectral methods.
 The code is open-source and developed by a team of researchers studying astrophysical, geophysical, and biological fluid dynamics.
 
-Dedalus is written primarily in Python and features an easy-to-use interface with symbolic vectorial equation entry.
-For example, to simulate incompressible hydrodynamics in a disk or ball, you can symbolically specify the PDEs, including boundary conditions and gauge constraints as well as the [tau modifications](https://dedalus-project.readthedocs.io/en/latest/pages/tau_method.html) needed to enforce them, as:
+Dedalus is written primarily in Python and features an easy-to-use interface with symbolic vectorial equation specification.
+For example, to simulate incompressible hydrodynamics in a ball, you can symbolically enter the equations, including [gauge constraints](https://dedalus-project.readthedocs.io/en/latest/pages/gauge_constraints.html) and [boundary conditions enforced with the tau method](https://dedalus-project.readthedocs.io/en/latest/pages/tau_method.html), as:
 
 ```python
-    problem.add_equation("div(u) + tau_c = 0")
-    problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_1) = - dot(u,grad(u))")
+    problem.add_equation("div(u) + tau_p = 0")
+    problem.add_equation("dt(u) - nu*lap(u) + grad(p) + lift(tau_u) = - dot(u,grad(u))")
     problem.add_equation("u(r=1) = 0")
     problem.add_equation("integ(p) = 0")
 ```
 
-Our numerical algorithms produce sparse and spectrally accurate discretiations of a wide variety of equations and domains, including Cartesian domains, disks, annuli, spheres, spherical shells, and balls, as seen in the following examples:
+Our numerical algorithms produce sparse and spectrally accurate discretizations of PDEs on simple domains, including Cartesian domains of any dimension, disks, annuli, spheres, spherical shells, and balls:
 
 <table style="background-color:#FFFFFF;">
   <tr>

docs/pages/changes_from_v2.rst

@@ -119,6 +119,15 @@ A best-of-both-worlds approach is to pass your scripts entire local namespace to
 This is achieved by passing the keyword ``namespace=locals()`` when instantiating problem objects.
 See the built in examples for illustrations of this approach to equation construction.
 
+Gauge constraints
+-----------------
+
+In Dedalus v2, gauge constraints (like the pressure gauge in incompressible hydrodynamics) were usually set by changing the equations for certain modes with the ``condition`` keyword when entering equations.
+In Dedalus v3, it's recommended to instead add spatially-constant gauge variables to the equations to introduce degrees of freedom that allow the gauge constraints to be directly imposed alongside the other equations.
+In most cases, the ``condition`` keyword can still be used if desired, but for technical/performance reasons it is no longer available in fully-Fourier problems.
+In any event, we find that the new approach (with gauge variables instead of conditions) makes for more readable equations.
+See the :doc:`gauge_constraints` page and the examples for more details.
+
 Tau terms
 ---------
 

docs/pages/gauge_constraints.rst

@@ -0,0 +1,59 @@
+Gauge Constraints
+*****************
+
+When you enter a system of PDEs in Dedalus, the left-hand side (LHS) of the equations is parsed into a sparse linear system.
+For the solver to succeed, this LHS matrix must be square and nonsingular.
+This means it must have a unique solution for *any* possible values on the right-hand side (RHS) of the equations, not just the particular RHS that you have entered.
+This also means that it must constrain *all* degrees of freedom of the variables, including gauge freedoms.
+
+For example, let's consider solving incompressible hydrodynamics in a fully periodic domain discretized with Fourier bases in each dimension.
+At first glance, we might think to build the problem variables and specify the equations simply as:
+
+.. code-block:: python
+
+    # Fields
+    p = dist.Field(name='p', bases=bases)
+    u = dist.VectorField(coords, name='u', bases=bases)
+
+    # Problem
+    problem = d3.IVP([p, u], namespace=locals())
+    problem.add_equation("div(u) = 0")
+    problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - dot(u,grad(u))")
+
+This formulation produces square matrices (same number of modes in the equations as the variables), but they are not all nonsingular.
+The problem is that the pressure gauge is undetermined -- any constant can be added to the pressure without affecting any of the equations.
+
+However, since the system is square, if one variable is underdetermined, then there must also be a degenerate constraint.
+In this case, the mean of the divergence equation is degenerate with the periodic boundary conditions imposed by the Fourier discretization.
+That is, for the mean :math:`\vec{k} = 0` Fourier mode, the divergence equation simply becomes ``"0 = 0"``.
+This seems consistent -- but again the system must be solvable for *any* RHS, not just the one we entered, and clearly we would have a problem if we entered ``"div(u) = 1"`` since for the mean mode this would become ``"0 = 1"``.
+
+To fix these problems, we need to add another equation that fixes the pressure gauge, and we need to remove the degenerate constraint.
+Algorithmically, we could simply replace the divergence constraint for the mean Fourier mode with a pressure constraint (as was done in Dedalus v2 using equation conditions).
+Another option is to expand the system by adding a spatially-constant variable (let's call it :math:`\tau_p`) to the divergence equation to absorb this degeneracy, and then impose the pressure gauge as a separate equation:
+
+.. code-block:: python
+
+    # Fields
+    p = dist.Field(name='p', bases=bases)
+    u = dist.VectorField(coords, name='u', bases=bases)
+    tau_p = dist.Field(name='tau_p')
+
+    # Problem
+    problem = d3.IVP([p, u, tau_p], namespace=locals())
+    problem.add_equation("div(u) + tau_p = 0")
+    problem.add_equation("dt(u) + grad(p) - nu*lap(u) = - dot(u,grad(u))")
+    problem.add_equation("integ(p) = 0")
+
+We've added a single additional degree of freedom to the variables and a single additional constraint, so this system is still square.
+It is now also nonsingular, since the mean pressure is fixed.
+The degeneracy of the mean divergence equation is also lifted as the tau variable can simply absorb/acquire the mean value of any possible RHS.
+The mean divergence of the velocity will always be zero, as required by the periodic discretization.
+
+Similar modifications work for other types of gauges and geometries.
+For example, for incompressible hydrodynamics in a bounded domain, we still need the above type of modification so that the integral of the divergence equation is compatible with the specified inflow boundary conditions.
+If the prescribed net inflow is nonzero, then the tau variable will acquire a corresponding nonzero value.
+From the modified equation, we can see that the velocity will then have a spatially uniform convergence equal to this tau value.
+Of course, for properly specified boundary conditions with no net inflow, the tau variable will be zero and the velocity will be divergence free.
+
+See the included :doc:`example scripts <tutorials>` for more examples of gauge modifications in various domains.

docs/pages/tau_method.rst

@@ -113,8 +113,7 @@ since
     \nabla \cdot G = \nabla^2 \vec{u} - \vec{\tau}_1(x) \partial_y P(y)
 
 Let's walk through setting up such a problem in Dedalus v3, assuming we're discretizing :math:`x` and :math:`y` with Fourier and Chebyshev bases, respectively.
-First, we need to create the necessary problem variable fields, including fields for the tau terms.
-We will also include a constant scalar tau term that will allow us to impose the pressure gauge.
+First, we need to create the necessary problem variable fields, including fields for the tau variables and a constant scalar tau for imposing the pressure gauge (see the :doc:`gauge_constraints` page):
 
 .. code-block:: python
 
@@ -137,8 +136,9 @@ Here we'll take :math:`P(y)` to be the highest mode in the Chebyshev-U basis, in
     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_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).
+We can then create a problem and enter the tau-modified PDEs, boundary condtions, and pressure gauge in vectorial form using these substitutions.
+Note that we will need to add the contant tau variable to the divergence equation as described in the :doc:`gauge_constraints` page.
+This allows us to impose the pressure gauge and removes the redundancy between the integral of the divergence equation and the integral of the inflow boundary conditions.
 
 .. code-block:: python
 
@@ -169,8 +169,8 @@ For instance, to enter the above equation set with homogeneous Dirichlet boundar
     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 disk and ball bases are not direct-product bases, so the tau terms can't actually be written just as the tau variable times a radial polynomial.
+Instead, for each horizontal mode (azimuthal mode :math:`m` in the disk and spherical harmonic :math:`\ell` in the ball), that mode of the tau variable 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.
 We've found that using tau polynomials from the original bases seems to give good results in the disk and ball:
 
@@ -202,7 +202,7 @@ To summarize, the main points regarding tau formulations are:
 3. For problems in Cartesian geometries, annuli, and spherical shells, we recommend a first-order-style implementation of the tau terms. Note that this only requires defining first-order substitutions that include tau terms, rather than increasing the problem size with first-order variables, as in Dedalus v2.
 4. For problems in the disk and ball, only a single tau term is needed for second-order elliptic/parabolic problems, and no first-order substitutions are necessary.
 
-See the included :doc:`example scripts <tutorials>` for more examples of working tau terms.
+See the included :doc:`example scripts <tutorials>` for more examples of tau modifications in various domains.
 
 
 .. .. math::

docs/pages/troubleshooting.rst

@@ -6,19 +6,8 @@ Singular matrix errors
 
 If you come across an error in the linear solver stating that a matrix/factor is singular, that means that the linear LHS portion of the PDE system is not uniquely solvable.
 This error indicates that some degrees of freedom of the solution are unconstrained and some of the specified equations are redundant (these are equivalent since the LHS matrices must be square).
-
 These errors are often due to imposing boundary conditions that are redundant for some set of modes and/or failing to constrain a gauge freedom in the solution.
-For instance, when simulating incompressible hydrodynamics, the pressure has a constant gauge freedom and the integral of the divergence constraint is redundant with the velocity boundary conditions.
-To fix these issues, a spatially constant variable can be added to the LHS of the divergence constraint, introducing a degree of freedom that allows for fixing the pressure gauge:
-
-.. code-block:: python
-
-    c = dist.Field() # Constant variable for fixing pressure gauge
-    ...
-    problem.add_equation("div(u) + c = 0") # Add c to break redundancy with boundary conditions
-    problem.add_equation("integ(p) = 0") # Pressure gauge
-
-See the :doc:`tau_method` page for more information on specifying boundary conditions.
+See the :doc:`gauge_constraints` and :doc:`tau_method` pages for more information on fixing these issues.
 
 Out of memory errors
 ====================

docs/pages/user_guide.rst

@@ -16,5 +16,6 @@ Specific how-to's:
 .. toctree::
     :maxdepth: 1
 
+    gauge_constraints
     tau_method
 

examples/ivp_2d_shear_flow/shear_flow.py

@@ -48,7 +48,7 @@
 p = dist.Field(name='p', bases=(xbasis,zbasis))
 s = dist.Field(name='s', bases=(xbasis,zbasis))
 u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis))
-tau_p = dist.Field(name='g')
+tau_p = dist.Field(name='tau_p')
 
 # Substitutions
 nu = 1 / Reynolds