fokker_planck_calculus

Module for functions used in calculating the integrals and doing the numerical differentiation for the implementation of the the full-F Fokker-Planck collision operator moment_kinetics.fokker_planck.

Parallelisation of the collision operator uses a special 'anyv' region type, see Collision operator and anyv region.

Struct to store the elemental nonlinear stiffness matrices used to express the finite-element weak form of the collision operator. The arrays are indexed so that the contraction in the assembly step is carried out over the fastest accessed indices, i.e., for YY0perp[i,j,k,iel], we contract over i and j to give data for the field position index k, all for the 1D element indexed by iel.

Struct to contain the integration weights for the boundary points in the (vpa,vperp) domain.

Struct of dummy arrays and precalculated coefficients for the Fokker-Planck collision operator when the Rosenbluth potentials are computed everywhere in (vpa,vperp) by direct integration. Used for testing.

Struct used for storing the integration weights for the boundary of the velocity space domain in (vpa,vperp) coordinates.

Struct of dummy arrays and precalculated coefficients for the finite-element weak-form Fokker-Planck collision operator.

Struct to store the boundary data for all of the Rosenbluth potentials required for the calculation.

Struct to contain data needed to create a sparse matrix.

Struct to store the (vpa,vperp) boundary data for an individual Rosenbluth potential.

Same as elliptic_solve!() above but no Dirichlet boundary conditions are imposed, because the function is only used where the lu_object_lhs is derived from a mass matrix. The source is made of two different terms with different weak matrices because of the form of the only algebraic equation that we consider.

Note: algebraic_solve!() run only in serial. They do not handle shared-memory parallelism themselves. The calling site must ensure that algebraic_solve!() is only called by one process in a shared-memory block.

Function to allocate an instance of vpa_vperp_boundary_data.

Function to allocate a boundary_integration_weights_struct.

Function to allocate at fokkerplanck_boundary_data_arrays_struct.

Function to allocate an instance of rosenbluth_potential_boundary_data.

Function to allocate an instance of sparse_matrix_constructor.

Function to assemble data in an instance of sparse_matrix_constructor. Instead of writing data.SS[icsc] = ss, as in assign_constructor_data!() we write data.SS[icsc] += ss.

Function to assemble the RHS of the kinetic equation due to the collision operator, in weak form. Once the array rhsvpavperp contains the assembled weak-form collision operator, a mass matrix solve still must be carried out to find the time derivative of the distribution function due to collisions. This function uses a purely parallel algorithm and may be tested by comparing to assemble_explicit_collision_operator_rhs_serial!(). The inner-most loop of the function is in assemble_explicit_collision_operator_rhs_parallel_inner_loop().

Function to assemble the RHS of the kinetic equation due to the collision operator, in weak form, when the distribution function appearing the derivatives is known analytically. The inner-most loop of the function is in assemble_explicit_collision_operator_rhs_parallel_analytical_inputs_inner_loop().

The inner-most loop of the parallel collision operator assembly. Called in assemble_explicit_collision_operator_rhs_parallel!().

Function to assemble the RHS of the kinetic equation due to the collision operator, in weak form. Once the array rhsvpavperp contains the assembled weak-form collision operator, a mass matrix solve still must be carried out to find the time derivative of the distribution function due to collisions. This function uses a purely serial algorithm for testing purposes.

Function to contruct the global sparse matrices used to solve the elliptic PDEs for the Rosenbluth potentials. Uses a dense matrix construction method. The matrices are 2D in the compound index ic which indexes the velocity space labelled by ivpa,ivperp. Dirichlet boundary conditions are imposed in the appropriate stiffness matrices by setting the boundary row to be the Kronecker delta (0 except where ivpa = ivpap and ivperp = ivperpp). Used for testing.

Function to contruct the global sparse matrices used to solve the elliptic PDEs for the Rosenbluth potentials. Uses a sparse matrix construction method. The matrices are 2D in the compound index ic which indexes the velocity space labelled by ivpa,ivperp. Dirichlet boundary conditions are imposed in the appropriate stiffness matrices by setting the boundary row to be the Kronecker delta (0 except where ivpa = ivpap and ivperp = ivperpp). See also assemble_matrix_operators_dirichlet_bc().

Function to assign data to an instance of sparse_matrix_constructor.

Function to assign precomputed (exact) data to an instance of vpa_vperp_boundary_data. Used in testing.

Function to allocated an instance of YY_collision_operator_arrays. Calls get_QQ_local!() from gauss_legendre. Definitions of these nonlinear stiffness matrices can be found in the docs for get_QQ_local!().

Function to carry out the direct integration of a formal definition of one of the Rosenbluth potentials, on the boundaries of the (vpa,vperp) domain, using the precomputed integration weights with dimension 4. The result is stored in an instance of vpa_vperp_boundary_data. Used in testing.

Function to carry out the direct integration of a formal definition of one of the Rosenbluth potentials, on the boundaries of the (vpa,vperp) domain, using the precomputed integration weights with dimension 3. The result is stored in an instance of vpa_vperp_boundary_data.

Function to carry out the integration of the revelant distribution functions to form the required coefficients for the full-F operator. We assume that the weights are precalculated. The function takes as arguments the arrays of coefficients (which we fill), the required distributions, the precomputed weights, the indicies of the `field' velocities, and the sizes of the primed vpa and vperp coordinates arrays.

Function to call direct integration function calculate_boundary_data!() and assign data to an instance of rosenbluth_potential_boundary_data, in place, without allocation.

Function to use the multipole expansion of the Rosenbluth potentials to calculate and assign boundary data to an instance of rosenbluth_potential_boundary_data, in place, without allocation. Use the exact results for the part of F that can be described with a Maxwellian, and the multipole expansion for the remainder.

Function to assign data to an instance of rosenbluth_potential_boundary_data, in place, without allocation. Used in testing.

Function to use the multipole expansion of the Rosenbluth potentials to calculate and assign boundary data to an instance of rosenbluth_potential_boundary_data, in place, without allocation.

Function to solve the appropriate elliptic PDEs to find the Rosenbluth potentials. First, we calculate the Rosenbluth potentials at the boundary with the direct integration method. Then, we use this data to solve the elliptic PDEs with the boundary data providing an accurate Dirichlet boundary condition on the maximum vpa and vperp of the domain. We use the sparse LU decomposition from the LinearAlgebra package to solve the PDE matrix equations.

Wrapper function to create a sparse matrix with an instance of sparse_matrix_constructor and sparse().

Elliptic solve function.

field: the solution
source: the source function on the RHS
boundary data: the known values of field at infinity
lu_object_lhs: the object for the differential operator that defines field
matrix_rhs: the weak matrix acting on the source vector
vpa, vperp: coordinate structs

Note: all variants of elliptic_solve!() run only in serial. They do not handle shared-memory parallelism themselves. The calling site must ensure that elliptic_solve!() is only called by one process in a shared-memory block.

Sets f(vpa,vperp) to a specied value f_bc at the boundaries in (vpa,vperp). f_bc is a 2D array of (vpa,vperp) where only boundary data is used. Used for testing.

Sets f(vpa,vperp) to a specied value f_bc at the boundaries in (vpa,vperp). f_bc is an instance of vpa_vperp_boundary_data.

Function to enforce boundary conditions on the collision operator result to be consistent with the boundary conditions imposed on the distribution function.

Unused function. Sets f(vpa,vperp) to zero at the boundaries in (vpa,vperp).

Function for getting the indices used to choose the integration quadrature.

get_global_compound_index(vpa,vperp,ielement_vpa,ielement_vperp,ivpa_local,ivperp_local)

For local (within the single element specified by ielement_vpa and ielement_vperp) indices ivpa_local and ivperp_local, get the global index in the 'linear-indexed' 2d space of size (vperp.n, vpa.n) (as returned by ic_func).

Function to get the local grid and integration weights assuming no divergences of the function on the 1D element. Gauss-Legendre quadrature is used for the entire element.

Function to get the local integration grid and quadrature weights to integrate a 1D element in the 2D representation of the velocity space distribution functions. This function assumes that there is a divergence at the point coord_val, and splits the grid and integration weights appropriately, using Gauss-Laguerre points near the divergence and Gauss-Legendre points away from the divergence.

ic_func(ivpa::mk_int,ivperp::mk_int,nvpa::mk_int)

Get the 'linear index' corresponding to ivpa and ivperp. Defined so that the linear index corresponds to the underlying layout in memory of a 2d array indexed by [ivpa,ivperp], i.e. for a 2d array f2d:

• size(f2d) == (vpa.n, vperp.n)
• For a reference to f2d that is reshaped to a vector (a 1d array) f1d = vec(f2d) than for any ivpa and ivperp it is true that f1d[ic_func(ivpa,ivperp)] == f2d[ivpa,ivperp].

Function that returns the sparse matrix index used to directly construct the nonzero entries of a 2D assembled sparse matrix.

Function to find the element in which x sits.

Function that precomputes the required integration weights only along the velocity space boundaries. Used as the default option as part of the strategy to compute the Rosenbluth potentials at the boundaries with direct integration and in the rest of (vpa,vperp) by solving elliptic PDEs.

Function that precomputes the required integration weights in the whole of (vpa,vperp) for the direct integration method of computing the Rosenbluth potentials.

Function to interpolate f(vpa,vperp) from one velocity grid to another, assuming that both grids are represented by (vpa,vperp) in normalised units, but have different normalisation factors defining the meaning of these grids in physical units. E.g.,

 vpai, vperpi = ci * vpa, ci * vperp
 vpae, vperpe = ce * vpa, ce * vperp

with ci = sqrt(Ti/mi), ce = sqrt(Te/mi)

scalefac = ci / ce is the ratio of the two reference speeds.

ivpa_func(ic::mk_int,nvpa::mk_int)

Get the vpa index ivpa that corresponds to a 'linear index' ic that spans a 2d velocity space.

Defined so that ivpa_func(inc_func(ivpa,ivperp,nvpa), nvpa) == ivpa.

See also ic_func, ivperp_func.

ivperp_func(ic::mk_int,nvpa::mk_int)

Get the vperp index ivperp that corresponds to a 'linear index' ic that spans a 2d velocity space.

Defined so that ivperp_func(inc_func(ivpa,ivperp,nvpa), nvpa) == ivperp.

See also ic_func, ivpa_func.

Base level function for computing the integration kernals for the Rosenbluth potential integration. Note the definitions of ellipe(m) ($E(m)$) and ellipk(m) ($K(m)$). https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.ellipe https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.ellipk

\[E(m) = \int^{\pi/2}_0 \sqrt{ 1 - m \sin^2(\theta)} d \theta\]

\[K(m) = \int^{\pi/2}_0 \frac{1}{\sqrt{ 1 - m \sin^2(\theta)}} d \theta\]

Function for computing the quadratures and carrying out the loop over the primed vpa coordinate in doing the numerical integration. Splits the integrand into three pieces – two which use Gauss-Legendre quadrature assuming no divergences in the integrand, and one which assumes a logarithmic divergence and uses a Gauss-Laguerre quadrature with an (exponential) change of variables to mitigate this divergence.

Function for computing the quadratures and carrying out the loop over the primed vpa coordinate in doing the numerical integration. Uses a Gauss-Legendre quadrature assuming no divergences in the integrand.

Function for computing the quadratures and carrying out the loop over the primed vperp coordinate in doing the numerical integration. Splits the integrand into three pieces – two which use Gauss-Legendre quadrature assuming no divergences in the integrand, and one which assumes a logarithmic divergence and uses a Gauss-Laguerre quadrature with an (exponential) change of variables to mitigate this divergence. This function calls loop_over_vpa_elements_no_divergences!() and loop_over_vpa_elements!() to carry out the primed vpa loop within the primed vperp loop.

The function loop_over_vperp_vpa_elements_no_divergences!() was used for debugging. By changing the source where loop_over_vperp_vpa_elements!() is called to instead call this function we can verify that the Gauss-Legendre quadrature is adequate for integrating a divergence-free integrand. This function should be kept until we understand the problems preventing machine-precision accurary in the pure integration method of computing the Rosenbluth potentials.

Function returns 1 for nelement > ielement >= 1, 0 for ielement = nelement.

Function returns 1 for nelement >= ielement > 1, 0 for ielement = 1.

Function returns 1 if igrid = ngrid or 0 if 1 =< igrid < ngrid.

Function returns 1 if igrid = 1 or 0 if 1 < igrid <= ngrid.

Function for getting the basic quadratures used for the numerical integration of the Lagrange polynomials and the integration kernals.

Function to compute the maximum error ${\rm MAX}|f_{\rm numerical}-f_{\rm exact}|$ for instances of vpa_vperp_boundary_data.

Function to compare two instances of rosenbluth_potential_boundary_data – one assumed to contain exact results, and the other numerically computed results – and compute the maximum value of the error. Calls test_boundary_data().