We seem to have an issue when adding a constant to a boundary-restricted field in problems with multidimensional bases. Here is an example inspired by a post to the users list:

import numpy as np
import dedalus.public as d3

coords = d3.PolarCoordinates('phi', 'r')
dist   = d3.Distributor(coords, dtype=np.float64)
disk   = d3.DiskBasis(coords, shape=(16, 16), radius=1, dtype=np.float64)
cb     = dist.Field(name='cb', bases=disk.edge)

(1+cb).evaluate()

Running this produces an error ValueError: non-broadcastable output operand with shape (1,1) doesn't match the broadcast shape (16,1). This is because the out field in AddFields has shape (1,1) even though its shape should be (16,1). Note that if the operator is cb+1, it code runs. The issue seems to be when building the bases for the out field. This occurs in Add’s _build_bases function. This loops over the coords in args[0].domain.bases_by_coord, and then adds up the bases of the two arguments coord by coord. The issue is that args[0] is 1, and its coords are the PolarCoordinates. However, cb's coords are phi. This means that the bases of the output are also None, which yields the shape of (1,1).

Is the issue that the 1 is not being converted to the proper field? Or is the behavior of _build_bases not expected?