CartesianGradient / Laplacian bug about dedalus HOT 4 CLOSED

Under further investigation, it appears that there is a term Laplacian(Ta) where Ta is a z-dependent NCC which is not being evaluated correctly. In commit ff3aea3 it is evaluated correctly, but in f263415 it is zero. Here's a simple code snippet which illustrates the problem:

import numpy as np
import dedalus.public as d3
coords = d3.CartesianCoordinates('x', 'z')
dist = d3.Distributor(coords, dtype=np.float64)
xbasis = d3.RealFourier(coords['x'], size=16, bounds=(0, 1), dealias=3/2)
zbasis = d3.ChebyshevT(coords['z'], size=16, bounds=(0, 1), dealias=3/2)
rhoa = dist.Field(name='rhoa', bases=(zbasis)) # isentropic density profile
Ta = np.exp(rhoa)
Lapa = d3.lap(Ta)
print(Lapa)

For ff3aea3 it returns Lap(exp(rhoa)) -- this is the desired behavior
For f263415 it returns 0

from dedalus.

Actually turns out the commit that causes the shift from Lap(exp(rhoa)) to 0 is 191879c, the commit before f263415.

from dedalus.

The issue with 191879c is this code here in the Laplacian _preprocess_args method:

elif not isinstance(coordsys, coords.DirectProduct) and operand.domain.get_basis(coordsys) is None:
    raise SkipDispatchException(output=0)

Since rhoa is a function of z only, rhoa.domain.get_basis(coordsys) returns None.

Presumably we want to be able to take the laplacian of variables which do not depend on all of the coordinates in a coordinatesystem?

from dedalus.

I verified that if I comment out those two lines in Laplacian _preprocess_args, the script acts as expected.

To summarize, it appears that 1218c62 introduced a bug in CartesianGradient which was fixed in f263415. However, the previous commit (191879c) introduced this "bug" into Laplacian.

We should carefully consider how Laplacian should act on fields which are constant in certain directions.

from dedalus.