Wall boundary condition experiment

[mrhardman]

Pull request for experimental feature for modifying the wall boundary condition.

[mrhardman]

Some figures demonstrating the potential usefulness of this feature for making sure that the kinetic-Chodura condition is satisfied without disturbing much the form or scaling of the potential with distance to the wall target (numerical dissipation in vpa only, not a full Fokker-Planck case).
wall_boundary_cutoff_scan.pdf

Generalising this cut-off to moment-kinetic or 2D cases is still left to be done. This feature may be useful for physical studies with the correct Fokker-Planck collision operator.

[johnomotani]

I'm wanting to try playing with this feature now, so I've merged in master. I thought it would be easier to merge if we store the boundary condition parameters within the coordinate structs, so there are fewer changes to function arguments - @mrhardman if you don't like that we can rethink. I've also added the feature to the moment-kinetic version of the sheath entrance boundary condition.

[mrhardman]

On balance I think I preferred my way. I think that putting the data in coordinate structs will force duplication of parameters, because every coordinate now will have this struct of data even if it is never used for anything (like in the initial conditions structs, where parameters that made sense for the spatial initialisation are also specified for the velocity initialisation, even if they are not used). In addition, if we just have the one parameter for now, making an extra name list that then gets put into coordinates struct seems like overkill when you could have just added an extra parameter to the coordinate inputs. I would rather store data according to theme (having the data in the boundary distributions struct, e.g.), and read the name list in from the boundary_conditions module. As the boundary distributions struct is already used in boundary_conditions.jl, I would hope that not too many extra arguments would need passing.

I can see why you want to go this way, and if you are taking ownership of the feature and strongly disagree with my points above I won't stop you : ). If we want to add similar features as this elsewhere, porting this to the boundary_conditions module might be something that I would want to do.

EDIT: I imagine examples of future parameter features would mostly be connected with how radial boundaries are treated, having a Maxwellian upstream and parameters to mimic the damping towards a vacuum downstream.

[mrhardman]

P.s. publication_inputs/xarray_post_processing/plot_wall.py from https://github.com/mabarnes/moment_kinetics/pull/228/files#diff-0768e749529d60b8847e50149f1e15df02852fc44d787e0b9f4cb0c51d587fe9 is the python script used to produce the .pdf plots of the epsz scan above. You might find the python script useful for making similar scans.

[johnomotani]

I see your point. When I started converting, I hadn't realised that the boundary_parameters were tucked away inside boundary_distributions, so there's not more to pass around. I don't have a strong preference, but just before I change it back - the main reason I see for preferring to put epsz in the coordinate is that we then keep the boundary conditions settings more or less in one place. Following that through, we should also put the bc setting in coordinate.boundary_parameters instead of as an individual field in coordinate. Would that consistency make you happier?

PS I've set it up so that extra parameters (only epsz so far!) are only created in the coordinate which actually uses them, so there are no redundant variables floating around.

[mrhardman]

The conservative option of keeping bc in coordinate would be my preference for now until we know what models we will need to address the radial boundary conditions, or more complex geometry with a private flux region as well as a LCFS, or even a domain that isn't neatly described by constant limits in coordinate values.

I say this partly because the coordinate inputs are still open to be refactored #241 into namelists, and partly because boundary condition options are potentially logically separate from the number and location of grid points to use (except in the finite difference method perhaps?). Perhaps we cannot reasonably take bc out of coordinate because of how periodic boundary conditions are dealt with, but I would prefer not to put anything into coordinate that makes it difficult (or inelegant) to make sandbox tests (e.g., as in https://github.com/mabarnes/moment_kinetics/blob/spatial_poissons_equation/test_scripts/PoissonSolve.jl, where I made a Poisson's equation solver in an infinite cylinder).

I am arguing for the conservative option of not changing coordinate, decide for yourself if my reasoning holds up : )

[johnomotani]

I'd think of boundary conditions as logically related to the coordinate, at least in the senses that: boundary conditions have to be specified for each coordinate, so it would make sense to do that in the input section for each coordinate; and in order to apply the boundary condition, you'd have to have the coordinate struct, at the very least to look up the irank and nrank for the distributed MPI.

I'm not against refactoring the coordinates input now to keep things simple. I think the only arguments that should be needed would be an OptionsDict and the name of the coordinate, so ad-hoc construction would look something like

foo_input = OptionsDict("foo" => OptionsDict("ngrid" => 5, "nelement" => 10))
define_coordinate(foo_input, "foo")

which would be simpler than it is at the moment because you could leave default values for options that you don't care about.

That refactor should probably be a separate PR, which we could then merge into this one (possibly via master). I'll make one and put it up as a proposal.

[mrhardman]

I suppose that I have to agree that information like "periodic" or "wall" or "none" is linked to the coordinate information in the way that the code in formulated, but I am not sure if very detailed information about distribution functions or numerical method choices (like epsz) should be included there, because this information is very specific to the solver in moment_kinetics and perhaps not very useful to keep in a general struct that can be used for other toy problems. Thanks for considering my point of view. I look forward to seeing the PR : ) Long live the namelists : )

[mrhardman]

@johnomotani Just left a comment in #241 pointing out a potential problem with the below approach.

foo_input = OptionsDict("foo" => OptionsDict("ngrid" => 5, "nelement" => 10))
define_coordinate(foo_input, "foo")

[johnomotani]

I think this feature is useful to have as an option. I also added some extensions to the post-processing in makie_post_processing that can show how far along in vpa we would have had to zero out f to satisfy the 'Chodura criterion'. Going to merge this once the tests pass.