Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow interpolation between grids when restarting #128

Merged
merged 21 commits into from
Nov 6, 2023

Conversation

johnomotani
Copy link
Collaborator

@johnomotani johnomotani commented Sep 11, 2023

Can be used to change resolution, and/or to switch between different moment-kinetic modes, when restarting.
Edit: Can also switch (in a basic way, see the docs) from 1D to 2D, or from 1V to 2V/3V.

This PR also includes some refactoring of the types used to dispatch to different discretization methods for differentiation and interpolation - this is a tidy-up and completion of the refactoring started in #57. Now there is an 'abstract type' discretization_info, defined in moment_kinetics_structs, which allows chebyshev_info (now a sub-type of discretization_info) to be moved back into chebyshev.jl, and replaces the use of Bool to indicate finite-difference methods with a finite_difference_info type (also a subtype of discretization_info).

TODO:

  • tests!
  • documentation

@johnomotani johnomotani added the enhancement New feature or request label Sep 11, 2023
@johnomotani johnomotani force-pushed the restart-interpolation branch 2 times, most recently from 4270457 to cd37b9f Compare September 13, 2023 09:19
@johnomotani
Copy link
Collaborator Author

[Copying this to the PR page because that will probably be easier to find in future!]

Does this commit allow for one to run a 1D1V input file to steady state (with r_ngrid = 1 and r_nelement = 1), and then restart it with r_ngrid > 1 and r_nelement > 1. This would be a great timesaver! 2D1V simulations might be more reliable if we start with a better initial condition -- which we might not be able to guess or write in terms of closed-form functions.

I haven't considered going from 1D to 2D. It should be possible, but will need a bit of work. The interpolation routines were written with velocity space interpolation in mind, so when they have to go past the end of the input grid they fill in values with an exp(-v^2) decay from the last value on the grid. I guess we'd want to add an option to control the extrapolation, which could be set for the r-dimension to something like 'use the last value on the grid'. I'm also not sure what exactly happens with a coordinate struct for a dimension with just one point. I remember that we do some sort of special-case for it, but that might break interpolation somehow. On the other hand, maybe we could take advantage of that to do a special-case interpolation method where if the dimension had only one point, then the 'interpolation' is just to copy the single value to the whole output array.

@mrhardman
Copy link
Collaborator

On the other hand, maybe we could take advantage of that to do a special-case interpolation method where if the dimension had only one point, then the 'interpolation' is just to copy the single value to the whole output array

This is what I had in mind. The special exception that you mention is in coordinates.jl -- grids with coord.n = 1 have their wgt[1] set to a constant and grid[1] set to 0.0. If we just plan to copy the data in the "radial dimension" from r.n = 1 to all the other indices, I would have thought we would have a fairly straight forward operation to do.

When reloading data for a restart, it is necessary to create
`coordinate` objects using the parallelisation of the new simulation,
not the parallelisation that was used to write the data.
Only the coordinate ranges are actually different between the
parallel_io=true and parallel_io=false branches, so can take a lot of
code outside the `if parallel_io` block.
This allows a simulation to restart from a run with different resolution
in any or all dimensions, and also to restart from a run with different
moment-kinetic settings (`evolve_density`, `evolve_upar` and
`evolve_ppar`).
This is needed to get the global_io_range correct, which is needed when
reloading.
This will allow them to be reused in other tests.
`discretization_info` is the abstract type inherited by all
implementations of a discretization (i.e. `chebyshev_info` and
`finite_difference_info`).

Also adds a `finite_difference_info` that is used to dispatch to the
finite difference methods, instead of using `Bool` for that.
Length-1 dimensions have to be treated a bit specially. Derivatives in
those directions should be zero. 'Interpolation' to a grid with more
than one point: for spatial dimensions, assume the variable is constant;
for velocity dimensions, use a Maxwellian whose peak value is the value
on the length-1 grid, and whose width is 1 (in units of the reference
speed) - this means that the integral of the variable over velocity
space, after accounting for factors of pi^0.5 or pi^1.5 in the
normalization of distribution functions, is the same (up to
discretization error) after 'interpolating' as it was before.
The special "chebyshev_pseudospectral_vperp" is now handled by having a
special case for "vperp" in the "chebyshev_pseudospectral" setup, so the
default for "verp_discretization" should just be
"chebyshev_pseudospectral".
The case when vpa and vz are not the same direction needs some special
handling (for the neutrals) to convert between 1V and 3V. This is not
implemented yet, so error in this case when bzeta!=0.
Have added tests for Krook collision operator and restart interpolation.
When not using parallel I/O, should not change the irank/nrank loaded
from the original output files.

Also check whether current distributed parallelisation is the same as
the original distributed parallelisation (changing this is only
supported when using parallel I/O).
@johnomotani johnomotani force-pushed the restart-interpolation branch from 20bac6b to 5b2f40e Compare October 31, 2023 15:45
@johnomotani johnomotani force-pushed the restart-interpolation branch from 5b2f40e to f13e3e0 Compare October 31, 2023 15:58
@johnomotani
Copy link
Collaborator Author

This branch now includes tests (https://github.com/mabarnes/moment_kinetics/blob/f13e3e03da956e775236d1a6d458eacb31dfcc9b/test/restart_interpolation_tests.jl) and handles interpolation from 1D to 2D and from 1V to 2V or 3V (48dcc80).

Should be ready to merge now, I hope.

@johnomotani johnomotani requested a review from mrhardman October 31, 2023 16:06
Base automatically changed from Krook-collisions to master November 1, 2023 12:37
Now that the `discretization_info` abstract type is defined, can define
the specific types for each implementation in the module along with the
implementation. This makes more sense and simplifies the dependencies
between different modules.
The method for setting up MPI-capable HDF5 has changed in recent
versions of HDF5.jl. Update the docs to reflect the new method.
This makes it easier to read and to refer to, as the amount of
information in the file has increased.
@johnomotani johnomotani merged commit f765bd7 into master Nov 6, 2023
8 of 21 checks passed
@johnomotani johnomotani deleted the restart-interpolation branch November 6, 2023 12:12
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants