Closed mrhardman closed 1 month ago
dfdvperp_and_conserving_experiment_plots.zip In this attachment I include plots which show the effect of turning on and off two two switches:
dFdvperp = 0
, use ad-hoc numerical conservationdFdvperp = 0
, use ad-hoc numerical conservationdFdvperp = 0
, do not use ad-hoc numerical conservation The results are (note that Sdot is not positive definite)
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).
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:
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.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.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
dFdvperp = 0
artificially for vperp advection and numerical dissipation, as this is currently the default boundary condition.vperp
boundary condition, perhaps by refactoring the "zero" option to be the current "zero-no-regularity" option.Two new test cases demonstrating the success of the vperp_bc="zero-no-regularity"
option:
nu_alphae
by the choice of input nuii
).
Tests showing that using vperp_bc="zero-no-regularity"
yields stable simulations with numerical dissipation
Tests of the laplacian_derivative!()
function in calculus_tests.jl
are added.
@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: https://github.com/mabarnes/moment_kinetics/commit/3cff3eb120b0633e1b41a5801c410341732afdce
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.
Tests using commit 8bff6a5658e8e3012fb8c0bef81c01c02760fb22:
@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 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.
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!
@johnomotani Happy for the merge to go ahead now providing the last commit passes the tests.
Pull request permitting simulations with additional collision model (numerical method) options. Based on #224, consider merging #224 first.
dFdvperp = 0
regularity condition. Setvperp_bc = "zero-no-regularity"
to access this feature.use_conserving_corrections = false
to access this feature.By default we set
vperp_bc = "zero"
anduse_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:
vperp_bc = "zero"
functionality to NOT imposedFdvperp = 0
. @johnomotani @mabarnes @LucasMontoya4 opinions?dFdvperp = 0
.dFdvperp = 0
.