CartesianGradient / Laplacian bug

[lecoanet]

There seems to be a bug in CartesianGradient which was found by Stephane Labrosse & detailed here:

https://groups.google.com/u/2/g/dedalus-users/c/UcOEEzv25TI

He has a compressible convection script that, when run for several thermal times, should reach thermal equilibrium with the same Nusselt number at top and bottom boundary, and with constant flux. This worked correctly up through commit ff3aea3, then the following commit 1218c62 which modified the CartesianGradient operator breaks the script and causes the simulation to generate NaNs within 10 timesteps. The later commit f263415 further modifies the CarteisanGradient operator -- at this commit, the script runs but does not generate correct results, e.g., the Nusselt number at top and bottom are different, and the total flux is not constant.

[lecoanet]

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

[lecoanet]

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

[lecoanet]

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?

[lecoanet]

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.