Implicit pseudo-timestep for electrons #266
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This can be used to make the grid spacing coarser at large wpa, which should help to relax CFL constraints in moment-kinetic runs.
Only works when the dimension is not distributed (i.e. when `coord.nrank==1`).
`t_params.implicit_coefficient_is_zero` only has length `n_rk_stages`, so need to check `istage < n_rk_stages` before getting `t_params.implicit_coefficient_is_zero[n_rk_stages+1]`.
If we just defined the residual for the electron distribution function solve to be 'dg/dt=0', then we would be asking the solver (roughly) to find g such that 'dg/dt<tol' for some tolerance. However dg/dt has some typical timescale, τ, so 'dg/dt∼g/τ'. It would be inconvenient to have to define the tolerances taking τ (normalised to the reference sound crossing time) into account, so instead estimate the relevant timescale as 'sqrt(me/mi)*z.L', i.e. as the electron thermal crossing time at reference parameters. We pass this as 'dt' to electron_kinetic_equation_euler_update!() so that it multiplies the residual.
`electron_kinetic_equation_euler_update!()` represents an explicit pseudo-timestep, so the -dt*(electron_ppar - electron_ppar_previous_ion_step)/ion_dt evaluation should be done with the 'old' pseudotimestep's electron_ppar, not the current estimate of the 'new' pseudotimestep's electron_ppar.
Initial implementation - currently uses fixed timestep.
The negative step in the line search was supposed to help make the iteration more robust by giving another option to make the residual decreases, but sometimes seems to stall convergence. May be better to just not use it.
Allows using `electron_backward_euler!()` for the first kinetic electron initialisation phase (where both electron distribution function and electron parallel pressure are evolved together) and using the `implicit_electron_ppar` option.
Prevents problems with convergence.
Can take care of incrementing `stage_counter[]` in `reset_nonlinear_per_stage_counters!()`, so it does not have to be done separately for each solver. Also rename to `reset_nonlinear_per_stage_counters!()` from `reset_nonlinear_per_stage_counters()` to indicate that this function modifies its argument.
For now, no parallelism in r supported - in future, need to add special functionality for that, similar to the 'anyv' regions, to support 2D simulations.
…_backward_euler!()
Actually only global within a block - if the dimension is parallelised with distributed MPI, it is not actually the global matrix.
…c-implicit-electron_ppar-loworder.toml Only really meant to include Krook collisions - not quite sure what 'nu_ei' does...
…ectron_backward_euler!() ...does not seem to actually reduce iteration counts significantly.
Should be marginally more efficient.
The handling of periodic bc only works when not using distributed-MPI, so is not very useful. It is now disabled by default, but can be enabled by passing a flag to `setup_gausslegendre_pseudospectral()`
For calculating preconditioners, it can be useful to have an explicitly calculated matrix ``` dense_second_deriv_matrix = inv(mass_matrix) * K_matrix ``` which includes the inverted mass-matrix already. Because of the matrix inverse, dense_second_deriv_matrix is a dense matrix, so (for efficiency) should not be used unless absolutely necessary.
...to ensure that the value of qpar and dqpar/dz is always consistent with the current distribution function, ppar, and vth
Also includes some fixes for the Jacobian matrix calculations, and extends them to handle a few more options.
'constant' boundary condition may still not be consistent with moment constraints, but does now use the 'unnormalised' speed to set the incoming Maxwellian distribution, and the incoming distribution is normalised by n/vth. Also fix normalisation - the 1/sqrt(pi) is now included in the normalisation of f (and in the integration routines) so does not need to be divided out here.
`implicit_electron_advance` tries to do a single non-linear solve for the steady state, and does not work (yet). `implicit_electron_ppar` uses pseudo-timestepping to find the steady state, with a backward-Euler pseudo-timestep, and does work (at least sometimes), so should be the default method.
Also fix up the debug checks for the kinetic electron run.
For debugging it is useful to see the errors, with a backtrace. This commit adds an optional flag to the makie_post_processing input that can switch off the error handling.
Should only be loaded when present in the restart file.
Allows different settings to be used for different electron sources, as intended, and gets rid of some warnings when multiple ion sources are used, but electron sources use defaults.
johnomotani
force-pushed
the
implicit-electrons
branch
from
8924b47 to
b8d42fa
Compare
johnomotani
force-pushed
the
implicit-electrons
branch
from
b8d42fa to
b767d22
Compare
Using Givens rotations (following 'Algorithm 2' in Zou (2023)) avoids the need for least-squares minimisations at each iteration of the GMRES linear solver.
MPI.bcast() can communicate (almost?) any type of object, but that means that the type of its result is not necessarily known before communication happens, leading to type instability. Therefore prefer to use other MPI.jl functions that are type-stable. Use in-place MPI operations in a few more places to avoid possibility of allocating extra buffers. Fix function wrapping in nonlinear_solvers to avoid type instability. Putting the function in a variable inside an if..elseif..else before wrapping it confuses the compiler, so instead need to do the 'wrapping' separately for each case. Fix way inner loop counter is used to avoid type instability Remove wrapper functions in nonlinear_solvers to avoid type instability
`Ref` is not a concrete type, so a struct defined with `Ref` members is not type-stable. The concrete type (relevant to our usage) is `Base.RefValue`.
johnomotani
force-pushed
the
implicit-electrons
branch
from
d496305 to
0aa3a46
Compare
Reduces need for synchronizations and reduces possibilities for bugs.
johnomotani
force-pushed
the
implicit-electrons
branch
from
0aa3a46 to
8eb96a9
Compare
Nonlinear solver (JFNK) seems to work on macOS now.
Better electron timestep size for kinetic electron test
johnomotani
force-pushed
the
implicit-electrons
branch
from
8eb96a9 to
e75d6ac
Compare
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Electrons use a pseudo-timestep loop to find steady state at each ion stage. Now the pseudo-timestep uses an implicit (backward-Euler) method, with the timestep size adjusted to be as long as possible while not taking an excessive number of iterations in the JFNK solver for each step.
Only tested so far for one case, see
examples/kinetic-electrons/README.md.Electron physics and solver methods still need to be documented - that will be done in a future PR.