Optimal FourierSeries implementation of Helmholtz operator

[astanziola]

The Helmholtz operator needs to evaluate the following differential operator internally

$$ \nabla^2u + \frac{1}{\rho}(\nabla \rho \cdot \nabla u) $$

with the caveat that the differential operators need to be modified to account for the PML

At the moment, this is done by modifying the equation as

$$ \hat \nabla^2u + \frac{1}{\rho}(\nabla \rho \cdot \nabla u), \qquad \hat \nabla^2 = \sum_{\xi \in {x,y,z}} \frac{1}{\gamma_{\xi}} \frac{\partial}{\partial \xi}\frac{1}{\gamma_{\xi}} \frac{\partial}{\partial \xi} $$

(see this).

However, in the case of FourierSeries fields, the partial derivatives are not implemented correctly for even derivatives: in particular, the Nyquist frequency is not handled correctly. See Algorithm 1 here for more info on this regard.

Any other kind of implementation seems to work worse. One could use the identity

$$ \hat \nabla^2u + \frac{1}{\rho}(\nabla \rho \cdot \nabla u) =\frac{1}{\rho}[\nabla \rho \cdot \nabla u] $$

and implement the RHS using Algorithm 4, which corresponds to the heterog_laplacian operator of jaxdf. This also seems to reduce accuracy, although the accuracy of heterog_laplacian has not been tested yet.

This calls for further investigation.

[astanziola]

Following the CBS paper, it may be enough to add the PML as a generic absorption term, without having to split the derivative operators and the Laplacian! This would save quite a bit of FFTs.

While testing the accuracy of this, it would be nice to start looking for both speed and error using basil