Closed mrhardman closed 6 months ago
Development in #149 since https://github.com/mabarnes/moment_kinetics/commit/c3b6261115ba671cf338e88bbb762de121f59c82 seems to have broken the long-time stability of sheath simulations with either numerical dissipation or fokker-planck collision operator. Examples to follow.
diss_test.zip For example, switching to numerical dissipation for a speedy simulation, at long times we have numerical spikes and a crash of the simulation due to negative temperatures. The same qualitative behaviour is seen for nuii = 1.0
time_evolving_0D2V_diss.zip An obvious bug is found if we run with only numerical dissipation in vperp and vpa. No sources are specified, but it looks like the solution is growing! @johnomotani Have I made a silly error in the input?
It appears that the "chebyshev_pseudospectral" option for vpa and vperp gives plausible behaviour but the "gausslegendre_pseudospectral" option gives anti-diffusive growth. @johnomotani What have we changed there recently?
This .gif is made by using the input above and changing the discretization options for vpa and vperp.
I suspect the problem is that the 'boundary terms' here https://github.com/mabarnes/moment_kinetics/pull/149/files#diff-6af4b87e0934f1c86e512a7fbaf16a705fcf73131def91f308817abbf9f14f24R1016-R1017 have the wrong sign - at least if I flip the signs of those then the evolution looks pretty much like the Chebyshev-discretization one. I haven't checked what those signs should be from the weak form discretization. If those signs are wrong, then I guess other boundary terms likely are as well?
I suspect the problem is that the 'boundary terms' here https://github.com/mabarnes/moment_kinetics/pull/149/files#diff-6af4b87e0934f1c86e512a7fbaf16a705fcf73131def91f308817abbf9f14f24R1016-R1017 have the wrong sign - at least if I flip the signs of those then the evolution looks pretty much like the Chebyshev-discretization one. I haven't checked what those signs should be from the weak form discretization. If those signs are wrong, then I guess other boundary terms likely are as well?
That sounds like a problem. Better tests for the second derivatives would be needed. Can you give a line number? The link doesn't resolve to specific file. If this is the issue, I don't think it affects the collision operator, which would explain why the Fokker-Planck relaxation test still passes. That means there is at least one other important bug.
EDIT: Thinking clearly, I did make tests of the boundary terms in https://github.com/mabarnes/moment_kinetics/blob/merge_fkpl_collisions/test_scripts/GaussLobattoLegendre_test.jl. This should mean that the extreme boundary terms are correct.
Improvements to the explicit boundary terms for the gausslegendre_pseudospectral option are here: https://github.com/mabarnes/moment_kinetics/pull/149/commits/d306695e713bbd35c5dadc5c3fd7ecc194de7c66. This seems to solve the problem with numerical dissipation in isolation, but it does not seem to explain why sheath simulations are broken.
Using https://github.com/mabarnes/moment_kinetics/pull/149/commits/6dfd0d9e310ab9de8fc9e385d9e3dc325aeda934 we still see problems in distributed memory simulations, even with a simple diffusion operator in a 1D1V simulation. The .gif here shows how the distribution function calculated incorrectly. The .zip contains an example of a local simulation that appears to perform as would be expected.
The currently parallelisation of C[F,F] requires us to use distributed memory for z for a simulation to run in a time that is affordable. The fact that the collision operator tests pass suggests that the remaining bug is in the distributed memory part of the code, which is not tested at check-in.
time_evolving_1D2V_new.zip Uploads including example input files and output .gif/.pdfs for a 1D2V sheath problem with C[F,F] and wall boundaries, and the corresponding Maxwellian relaxation problem (0D2V). These simulations suggest that commit c3b6261115ba671cf338e88bbb762de121f59c82 behaves acceptably for 1) small enough time steps and 2) simulations where the stagnation point z=0 are eliminated with the choice of grids.
_Originally posted by @mrhardman in https://github.com/mabarnes/moment_kinetics/issues/149#issuecomment-1825623348_