Fokker Planck optional features

[mrhardman]

Pull request permitting simulations with additional collision model (numerical method) options. Based on #224, consider merging #224 first.

• Collisions that avoid enforcing the dFdvperp = 0 regularity condition. Set vperp_bc = "zero-no-regularity" to access this feature.
• Collisions that time-evolve without imposing the ad-hoc correction terms. Set use_conserving_corrections = false to access this feature.

By default we set vperp_bc = "zero" and use_conserving_corrections = true.

A check-in test is provided for vperp_bc = "zero-no-regularity". Additional experience is required to determine if this option is reliable for general simulations.

Decisions/Actions required:

• Whether or not to refactor vperp_bc = "zero" functionality to NOT impose dFdvperp = 0. @johnomotani @mabarnes @LucasMontoya4 opinions?
• Test vperp numerical diffusion with/without dFdvperp = 0.
• Test vperp advection without imposition of dFdvperp = 0.

[mrhardman]

dfdvperp_and_conserving_experiment_plots.zip
In this attachment I include plots which show the effect of turning on and off two two switches:

• fokker-planck-relaxation-beam-init1: impose dFdvperp = 0, use ad-hoc numerical conservation
• fokker-planck-relaxation-dfdvperp-test1: do not impose dFdvperp = 0, use ad-hoc numerical conservation
• fokker-planck-relaxation-dfdvperp-test1-no-conserve: do not impose dFdvperp = 0, do not use ad-hoc numerical conservation

The results are (note that Sdot is not positive definite)

• fokker-planck-relaxation-beam-init1: Sdot tends to zero, moments are conserved to ~ 1.0e-6 over 20 collision times
• fokker-planck-relaxation-dfdvperp-test1: Sdot tends to zero (but oscillates around zero, going negative), moments are conserved to ~ 1.0e-12 over 20 collision times
• fokker-planck-relaxation-dfdvperp-test1-no-conserve: Sdot goes negative, but approaches zero from below. Moments are not conserved, which changes of order unity by 20 collision times.

Overall this justifies the use of ad-hoc conservation terms for long-term stability, and raises the question of whether we need to use dFdvperp = 0 as an imposed boundary condition when the collision operator is implemented with FEM and has zero flux at the origin. dFdvperp = 0 was introduced to stabilise a simulation with a numerical diffusion which had a flux at vperp = 0 (i.e., it was d^2 F/dvperp^2).

[mrhardman]

Further tests with the cross-collision operator looking at the effect of imposing dFdvperp =0. In these tests no sources or sinks are used. There are three test cases here:

• sd_test_impose_dfdvperp_zero_large_dt.zip: Here we use numerical conserving terms, and we impose dFdvperp =0. We use a timestep that is too large, and the system goes numerically unstable around t = 14. Note the strange behaviour of the thermal speed at this time, and that the system seems to recover a physical-looking solution, even though we would rather find out more easily that this run is broken when we look at the final time step.
• sd_tests_large_dt.zip: Here we repeat the test above, without imposing dFdvperp =0 artificially in the time stepper. We still use numerical conserving terms, and the system still goes numerically unstable, but this time the system does not "recover", and we can easily tell that the solution is broken.
• sd_df_tests_small_dt.zip: Here we reduce the time step size, whilst we still use the "zero-no-regularity" boundary condition, meaning that dFdvperp =0 is not imposed. Now the simulation runs without numerical instability and the thermal speed converges to a fixed value (determined by the stationary background distributions that the evolved species collides with).

These results suggest that

• We should determine whether we ever need to impose dFdvperp = 0 artificially for vperp advection and numerical dissipation, as this is currently the default boundary condition.
• We should consider changing the default vperp boundary condition, perhaps by refactoring the "zero" option to be the current "zero-no-regularity" option.

[mrhardman]

Two new test cases demonstrating the success of the vperp_bc="zero-no-regularity" option:

• sd_df_source_sink_self_collisions_false.zip: A test set up with sufficient scale separation and losses at small v to find the canonical isotropic slowing down distribution function, with comparison against the analytical result. Self-collisions are neglected. The simulation runs for around 20 normalised times (note that the time scale is normalised to nu_alphae by the choice of input nuii).
• sd_df_source_sink_self_collisions_true.zip: The same test as above, with self collisions. Note that the slowing-down analytical solution does not match so well for small velocities, in accordance with expectations.

[mrhardman]

Tests showing that using vperp_bc="zero-no-regularity" yields stable simulations with numerical dissipation

• numerical_dissipation_cheb.zip: Test with Chebyshev pseudo spectral representation of Laplacian operator.
• numerical_dissipation_gausslegendre.zip: Test with Gauss Legendre finite element representation of Laplacian operator.

Tests of the laplacian_derivative!() function in calculus_tests.jl are added.

[johnomotani]

Whether or not to refactor vperp_bc = "zero" functionality to NOT impose dFdvperp = 0.

I think this sounds like a good idea - switch so that the default vperp_bc = "zero" does not impose dFdverp = 0 and the option that does impose it has a different name ("zero-impose-regularity" maybe? I'd leave it up to @mrhardman to pick a good name).

[mrhardman]

@johnomotani I think that the dirichlet_bc feature in the Gauss Legendre derivatives is breaking my tests. Are you using these arrays to solve 2nd order ODEs? If not, I would remove this and impose boundary conditions externally to the derivatives.

Assumed problem commit: 3cff3eb

EDIT: by making the dirichlet_bc option leave the lower endpoint of vperp alone, I can make my tests pass. However, I still would like some clarification on how these matrices are being used, as the boundary conditions should only be applied if K and L are being inverted. I don't expect this is the default. In my opinion, new matrices for inverting would be a better option.

[mrhardman]

Tests using commit 8bff6a5:

• advection_test_with_collisions_no_regularity.zip: A test showing that we do not need to impose dFdvperp = 0 for collisions to run with apparently physical behaviour (until size of vpa vperp domain is too small) with dvperp/dt advection terms.
• advection_test_no_collisions_no_regularity.zip: A test showing that we do not need to impose dFdvperp = 0 for apparently physical behaviour with dvperp/dt advection (until size of vpa vperp domain is too small).

[johnomotani]

@mrhardman I think I must've added that 'feature' to make matrices including the boundary condition that I could invert for preconditioning some implicit solves. I think I'll end up adding together several of these matrices, and modifying the result to impose Dirichlet bc anway, so I think the 'Dirichlet-bc-in-the-matrix' can be removed. If it's useful I think you're right and having separate matrices would be the right thing to do.

[mrhardman]

@mrhardman I think I must've added that 'feature' to make matrices including the boundary condition that I could invert for preconditioning some implicit solves. I think I'll end up adding together several of these matrices, and modifying the result to impose Dirichlet bc anway, so I think the 'Dirichlet-bc-in-the-matrix' can be removed. If it's useful I think you're right and having separate matrices would be the right thing to do.

How would you like to proceed? Should I revert the commits that made this feature, or do you want to deal with this in a new PR where you look at the implicit solves?

EDIT: Or I could just make a new commit changing the default, and let you add new matrices at a later date.

[johnomotani]

Or I could just make a new commit changing the default, and let you add new matrices at a later date.

Yeah, that sounds good!

[mrhardman]

@johnomotani Happy for the merge to go ahead now providing the last commit passes the tests.