Does CyclicLaplace actually work with periodicX=true?

[johnomotani]

I just tried running the Hasegawa-Wakatani example with periodicX = true (also using hw:modified=false, but I don't think that's relevant here). It did run, but very badly until I also set the INVERT_ZERO_DC flag. That makes me suspect that there is an issue with the DC component when periodicX = true. I would have thought it needed some special handling (e.g. using Dirichlet boundary conditions for the DC component, so that the kz=0, kx=0 piece is fixed), but don't see any in src/invert/laplace/impls/cyclic/cyclic_laplace.cxx or include/cyclic_reduction.hxx. There is some handling of 'periodic' that is the same for all kz's. Does anyone know whether Laplace solves with periodicX = true ought to work?

[bendudson]

@johnomotani I think it works in the sense that the code does calculate what it's designed to. Unfortunately the design has a problem: If both X and Z are periodic then the (0,0) mode is unconstrained, as it is when Neumann boundary conditions are applied to both boundaries, so the matrix being inverted for the kz=0 component is singular. It seems to have a go at a solution, but runs poorly because I guess the solution is very numerically unstable.

Setting INVERT_ZERO_DC removes the troublesome (0,0) mode, but unfortunately also removes all (kx, 0) modes too, which are the interesting zonal modes.

To fix this we would need to remove the null-space (0,0) mode. I think we probably want to remove the average phi, rather than pin an arbitrary point to zero, but doing that in the matrix inversion breaks the tridiagonal structure. For the kz = 0 mode we can probably fix an arbitrary X location to zero in the tridiagonal inversion, calculate phi, and then remove the X average...?

[johnomotani]

@bendudson sounds reasonable to me! :+1: @spamela was asking me about Hasegawa-Wakatani, so I'll leave it up to the interested people whether this is worth pushing on. Should be a fairly simple fix though.

[spamela]

Hey all, Yes, there is definitely interest in having properly periodic boundary conditions for the HW model. Let me know!

[johnomotani]

So if someone wants to work on this (i.e. make a branch off next and submit a PR), I think the main thing is that there needs to be a special case for localmesh->periodicX == true here https://github.com/boutproject/BOUT-dev/blob/d34cd444b6f63e6c9c2dc6fe543b07c00b1a04b0/src/invert/laplace/invert_laplace.cxx#L477 which handles the kz == 0 mode by setting matrix coefficients to force some point in the result equal to the rhs. Then after the solve here https://github.com/boutproject/BOUT-dev/blob/d34cd444b6f63e6c9c2dc6fe543b07c00b1a04b0/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx#L178 and here https://github.com/boutproject/BOUT-dev/blob/d34cd444b6f63e6c9c2dc6fe543b07c00b1a04b0/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx#L247 before the FFT back to real space do as @bendudson suggested and subtract the x-average from the kz == 0 mode (not sure we actually have a function to do that for a single Fourier mode - likely to need an MPI_Allreduce over the x-direction-communicator I guess). Finally there should be a test (part of or similar to https://github.com/boutproject/BOUT-dev/tree/next/tests/integrated/test-laplace) to verify that Laplace inversion works correctly in the periodicX case.

[mredenti]

I just tried running the Hasegawa-Wakatani example with periodicX = true (also using hw:modified=false, but I don't think that's relevant here). It did run, but very badly until I also set the INVERT_ZERO_DC flag. That makes me suspect that there is an issue with the DC component when periodicX = true. I would have thought it needed some special handling (e.g. using Dirichlet boundary conditions for the DC component, so that the kz=0, kx=0 piece is fixed), but don't see any in src/invert/laplace/impls/cyclic/cyclic_laplace.cxx or include/cyclic_reduction.hxx. There is some handling of 'periodic' that is the same for all kz's. Does anyone know whether Laplace solves with periodicX = true ought to work?

In the meantime, are there any other methods I can use to solve the Laplace inversion problem with periodicX = true?

[bendudson]

Hi @mir1995 I think that the options are:

• setting the INVERT_ZERO_DC flag, but understanding that the zonal flows will be missing. This may be important, depending on the regime.
• Not using a periodic domain: Just set periodicX=false, and use a larger box so that the region you're analysing is far from any boundary

Unfortunately other solvers, such as the PETSc implementation, would also need changes to work with periodic X boundaries. Sorry, periodic X is not an option that is used very often.

[mredenti]

Hi @mir1995 I think that the options are:

• setting the INVERT_ZERO_DC flag, but understanding that the zonal flows will be missing. This may be important, depending on the regime.
• Not using a periodic domain: Just set periodicX=false, and use a larger box so that the region you're analysing is far from any boundary

Unfortunately other solvers, such as the PETSc implementation, would also need changes to work with periodic X boundaries. Sorry, periodic X is not an option that is used very often.

Hi @bendudson, thank you for your reply. I would have thought that any finite difference scheme would have worked instead. Does the matrix A generated by the finite difference + boundary conditions in the system 'A\phi = \vort' become singular?

[bendudson]

Hi @mir1995 Yes that's the problem: In a doubly-periodic box an arbitrary constant offset can be added to the potential phi in Del^2 phi = vort, so the matrix is singular with a null space.

[mredenti]

Hi @mir1995 I think that the options are:

• setting the INVERT_ZERO_DC flag, but understanding that the zonal flows will be missing. This may be important, depending on the regime.

I guess one could solve separately for the zonal components and then add it to the solution? As a temporary fix that is.

[bendudson]

@mir1995 Difficult... the zonal components evolve along with the turbulence, so are not something that can be calculated once and added together: The zonal flows take energy from the turbulence, and interact with it in complex ways with a predator-prey kind of behaviour. I think fixing the solver as John wrote above should be fairly quick, and will try to have a go if I get time this week,

[mredenti]

@mir1995 Difficult... the zonal components evolve along with the turbulence, so are not something that can be calculated once and added together: The zonal flows take energy from the turbulence, and interact with it in complex ways with a predator-prey kind of behaviour. I think fixing the solver as John wrote above should be fairly quick, and will try to have a go if I get time this week,

@bendudson By averaging the modified HW equations in the z-direction one should obtain an equation for the zonal components which off course would need to be added to the physics model in order to be solved.

[bendudson]

@mir1995 Unfortunately I think solving that equation would require a modified 1D periodic solver, so doesn't get around the need to modify the cyclic inversion code.

[mredenti]

@mir1995 Unfortunately I think solving that equation would require a modified 1D periodic solver, so doesn't get around the need to modify the cyclic inversion code.

@bendudson ah! off course. Thank you for your reply.

[bendudson]

Hi @mir1995 @johnomotani @spamela I think this periodic X domain is now fixed in PR #2557 It's missing a test to check it though...