Skip to content

Implicit, Fokker-Planck, self collisions for ion species - #353

Closed
mrhardman wants to merge 169 commits into
masterfrom
implicit-fp-collisions-precon-debug
Closed

Implicit, Fokker-Planck, self collisions for ion species#353
mrhardman wants to merge 169 commits into
masterfrom
implicit-fp-collisions-precon-debug

Conversation

@mrhardman

@mrhardman mrhardman commented Apr 23, 2025

Copy link
Copy Markdown
Collaborator

The aim of this PR is to provide the ability to advance the equation

$$\frac{\partial F_i}{\partial t} = C[F_i,F_i] \quad (1)$$

with the implicit, backward-Euler timestepping method, i.e., solve

$$F^{n+1}_i - F^n_i = \Delta t C[F^{n+1}_i, F^{n+1}_i] \quad (2)$$

using the newton_solve!() functionality provided by @johnomotani. This feature is integrated into the main time advance functions of moment_kinetics so that the implicit implementation of the collision operator can be used with one of the IMEX time advance schemes also provided by @johnomotani.

An example input file that uses an IMEX advance to solve the equation (1) above is provided in examples/fokker-planck/fokker-planck-relaxation-implicit.toml. Note the addition of the input section [fokker_planck_collisions_nonlinear_solver] which provides the flags for newton_solve!(), and note the presence of a new @enum for selecting the type of ion advance, e.g.,

    [fokker_planck_collisions]
    use_fokker_planck = true
    nuii = 1.0
    frequency_option = "manual"

    [fokker_planck_collisions_nonlinear_solver]
    atol = 1.0e-10
    rtol = 1.0e-10
    
    [timestepping]
    nstep = 10
    dt = 1.0e-3
    nwrite = 1
    nwrite_dfns = 1
    kinetic_ion_solver = "implicit_ion_fp_collisions"
    type = "PareschiRusso3(4,3,3)"

A separate set of inputs is provided for the ion collision non-linear solve as for the electron non-linear solve in case different default inputs are necessary. At this time I have not tested implicit ion collisions alongside implicit kinetic electrons, but it is plausible that different defaults could be convenient for fast simulations. The list of possible kinetic_ion_solver options are given here https://github.com/mabarnes/moment_kinetics/blob/11e4f651bce559643aa2f179efb630fdf1c3cec2/moment_kinetics/src/input_structs.jl#L51-L65.

A series of test scripts have been developed in the building of this feature. In test_scripts/ImplicitCollisionsTest.jl we have the test scripts test_implicit_collisions() https://github.com/mabarnes/moment_kinetics/blob/11e4f651bce559643aa2f179efb630fdf1c3cec2/test_scripts/ImplicitCollisionsTest.jl#L63-L75 and test_implicit_collisions_wrapper() https://github.com/mabarnes/moment_kinetics/blob/11e4f651bce559643aa2f179efb630fdf1c3cec2/test_scripts/ImplicitCollisionsTest.jl#L217-L230, which can be used to solve equation one with backward-Euler and the initial condition

$$F_i(t=0) \propto exp\left(-\frac{(v_\| - v_{\| 0})^2 + (v_\perp - v_{\perp 0})^2}{v_{th0}^2}\right) \quad (3)$$

An example function call to step with one timestep of backward-Euler would be

julia> using Revise; using moment_kinetics
julia> include("test_scripts/ImplicitCollisionsTest.jl")
julia> test_implicit_collisions_wrapper(test_particle_preconditioner=true,test_numerical_conserving_terms=true,
    vth0=0.5,vperp0=1.0,vpa0=1.0, nelement_vpa=32,nelement_vperp=16,Lvpa=8.0,Lvperp=4.0, bc_vpa="none", bc_vperp="none",
     ntime=1, delta_t = 1.0, ngrid=5)

I have also extended the coverage of the automatic test scripts moment_kinetics/test/fokker_planck_tests.jl and moment_kinetics/test/fokker_planck_time_evolution_tests.jl. The changes in moment_kinetics/test/fokker_planck_tests.jl test

  • the preconditioner based on the backward-Euler solution of the equation
$$\frac{\partial F_i}{\partial t} = C[F_i,F_0] \quad (1)$$

with $F_0$ a specified distribution function, usually chosen to be $F_i$ at the last timestep;

  • and test the lowest level nonlinear backward Euler solve of equation (1) via equation (2).

The changes in moment_kinetics/test/fokker_planck_time_evolution_tests.jl

julia> using Revise; using moment_kinetics
julia> include("moment_kinetics/test/fokker_planck_time_evolution_tests.jl")
julia> FokkerPlanckTimeEvolutionTests.print_output_data_for_test_update("runs/path_to_your_run/path_to_your_run", write_grid=true, write_pdf=true)

In summary, this PR

  • adds the option for implicit backward-Euler advance with nonlinear ion self collisions,
  • provides a preconditioner that is developed from the linearised test-particle collision operator that might be used to make cheaper implicit model operators,
  • gives the option to use near-exact numerical error corrections to preserve the moments of $F_i$ (the default behaviour),
  • extends newton_solve!() to support iterations in the anyv region layout,
  • tests the feature with automatic tests,
  • and updates the tests to use the initial condition (3) to give a clear comparison in tests to published examples

Before merging, I still need to add

  • An error exception to catch use of vperp.bc="zero-impose-regularity", which is not supported.
  • Revert changes to set_defaults_and_check_values!() and make [fokker_planck_collisions.nonlinear_solver] a standalone namelist which is read directly from the input TOML in setup_fp_nl_solve(). The function setup_fp_nl_solve() should contain any Fokker-Planck specific defaults for the nonlinear solver. Remove the nonlinear_solver member from fkpl_collisions_input.
  • Include nl_solver_params.ion_fp_collisions in suitable timestep checks.
  • Set defaults to use preconditioner (currently preconditioner is not used apart from in test scripts) and supply sensible default way of choosing to update it (note that this requires some method of storing the preconditioner at all spatial points).

Before merging this PR I would request feedback on

  • structure of inputs to select the kinetic ion time advance options and the non-linear solver options
  • suitability of the tests
  • whether or not sufficient documentation is provided
  • whether checks with simulations with vpa, vperp, z advection are required
  • anything else that comes up in discussion.

mrhardman and others added 30 commits February 21, 2025 14:26
…), this removes random NaN/Inf errors that seem to have to do with these arrays not being assigned zero values.
…lthough still see Inf when using this option).
…, vpa) for testing implicit collision advance.
…d_dot for Newton-Krylov method. Now running test_scripts/ImplicitCollisionsTest.jl shows the same result (up to 1e-5) for runs with 1 and 2 cores.
…ne when using test_scripts/ImplicitCollisionsTest.jl non-interactively.
…ed when changing regions at line 285 of test_scripts/ImplicitCollisionsTest.jl. The error seems not to be fixed by synchronizing the anyv block within newton_solve!. The affected array appears to be Fnew, but obtaining a more helpful stacktrace is hindered because the anyv region synchronizations are the only ones permitted within newton_solve! in this setup, but only the serial_region macro appears to trigger the error.
…, test_implicit_collisions(ngrid=5,nelement_vpa=32,nelement_vperp=16,test_particle_preconditioner=true,test_linearised_advance=true,ntime=10,delta_t=5.0,Lvpa=8.0, Lvperp=4.0, vperp0=0.0).
…fied in test particle preconditioner matrix. A further refactor to introduce a single calculate_rosenbluth_potential!() wrapper may be beneficial.
…check that preconditioner matrix does indeed force F -> F_Maxwellian in steady state.
…, to run test shown in figure 8 of preprint, use test_implicit_collisions(ngrid=5,nelement_vpa=32,nelement_vperp=16, test_linearised_advance=false, test_particle_preconditioner=true, use_Maxwellian_Rosenbluth_coefficients_in_preconditioner=false, test_numerical_conserving_terms=true, delta_t=1.0, ntime=200, Lvpa=8.0, Lvperp=4.0, vth0=0.5). Note the order of magnitude speedup switching between test_particle_preconditioner=false and test_particle_preconditioner=true.
CI tests on Github Actions are failing due to segfaults in NCDatasets
calls. The parallel tests do not seem to be affected, don't know why -
maybe because they compile a system image??

Hopefully the upstream packages will sort this out eventually. To keep
an eye on when this happens, have kept an 'examples' test job that
includes NCDatasets, which is expected to fail at the moment. When that
job passes again, we should revert this commit.
Although other tests segfault with NCDatasets, HDF5 fails to precompile
when using OpenMPI_jll unless NCDatasets is also installed.
mrhardman added 2 commits May 9, 2025 13:36
…this is the only remaining test failure."

This reverts commit 1f26367.
…the pre-(#322)-merge state where the run time diagnostic for qpar used only the 1V definition and did not include the contribution from `vperp^2` to the energy transported by the flux. To see the corresponding change in the pre-merge code, see b987048 and #365.
@mrhardman

Copy link
Copy Markdown
Collaborator Author

I have not included implicit FP collisions in the timestep size checks
https://github.com/mabarnes/moment_kinetics/blob/d0b79f0687f81a21da715258df44ef2b2a129e67/moment_kinetics/src/time_advance.jl#L3422-L3467

. @johnomotani Should I do so by extending the lists of solver parameters?

Yes, I think so.

OK. I will address this after merging #363 into this PR.

…ster

Implicit FP collisions: merge master, changing normalisations to new conventions.
@mrhardman

Copy link
Copy Markdown
Collaborator Author

@johnomotani Although 167ce15 passed all tests, the tests are now failing post merge of #363, and I see still over a thousand files with changes with respect to master. I think the former could be explained by mistakes in automatic merge, although I have not checked what the error is yet. The latter problem seems more serious, since this PR should now include all the normalisation changes etc from #322. Do you see what has gone wrong?

johnomotani and others added 14 commits June 16, 2025 17:59
…ster

Implicit fp collisions merge master, further attempt.
Aiming to make performance tests using this script more representative:
* Parallelise array update
* Make some diagnostics optional
* Start 'run' timing immediately before time-advance loop to exclude all
  initialisation
* Lift struct lookups out of inner loops.
* Put @inbounds call in mass_matrix_interior_loop()
* Reduce number of multiplications in
  assemble_explicit_collision_operator_rhs_parallel_inner_loop()
It should do this as the 'anyv' subblocks are subblocks of the
shared-memory block that is synchronized by _block_synchronize(), so
every subblock is synchronized. Done by calling
_anyv_subblock_synchronize() - although this is mildly inefficient
(MPI.Barrier() gets called twice), it is only used for debugging so this
should be OK.
Use `comm_anyv_subblock[]` instead of the incorrect `comm_block[]`.
…ebug-parallelise-precon-construction

Parallelise Fokker-Planck preconditioner matrix construction
@johnomotani

Copy link
Copy Markdown
Collaborator

Replaced by #383.

@johnomotani
johnomotani deleted the implicit-fp-collisions-precon-debug branch July 13, 2025 15:51
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.

Implicit Fokker-Planck collision advance with JFNK solver

2 participants