From 8d27207f9dc21204799c3a4584166b43e5dd6d39 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Thu, 24 Oct 2024 15:36:22 +0100 Subject: [PATCH 1/6] Initial improvement of autodocs for fokker_planck.jl --- moment_kinetics/src/fokker_planck.jl | 84 +++++++++++++++++++--------- 1 file changed, 59 insertions(+), 25 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index fb5ac876b..521dca83d 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -1,5 +1,5 @@ """ -module for including the Full-F Fokker-Planck Collision Operator +Module for including the Full-F Fokker-Planck Collision Operator. The functions in this module are split into two groups. @@ -77,10 +77,10 @@ using ..reference_parameters: setup_reference_parameters Function for reading Fokker Planck collision operator input parameters. Structure the namelist as follows. -[fokker_planck_collisions] -use_fokker_planck = true -nuii = 1.0 -frequency_option = "manual" + [fokker_planck_collisions] + use_fokker_planck = true + nuii = 1.0 + frequency_option = "manual" """ function setup_fkpl_collisions_input(toml_input::Dict) reference_params = setup_reference_parameters(toml_input) @@ -152,7 +152,7 @@ end ######################################################## """ -function that initialises the arrays needed for Fokker Planck collisions +Function that initialises the arrays needed for Fokker Planck collisions using numerical integration to compute the Rosenbluth potentials only at the boundary and using an elliptic solve to obtain the potentials in the rest of the velocity space domain. @@ -238,7 +238,8 @@ end """ Function for advancing with the explicit, weak-form, self-collision operator using the existing method for computing the Rosenbluth potentials, with -the addition of cross-species collisions against fixed Maxwellian distribution functions +the addition of cross-species collisions against fixed Maxwellian distribution functions +where the Rosenbluth potentials are specified using analytical results. """ @timeit global_timer explicit_fp_collisions_weak_form_Maxwellian_cross_species!( pdf_out, pdf_in, dSdt, composition, collisions, dt, @@ -305,7 +306,7 @@ end """ -Function for advancing with the explicit, weak-form, self-collision operator +Function for advancing with the explicit, weak-form, self-collision operator. """ @timeit global_timer explicit_fokker_planck_collisions_weak_form!( pdf_out, pdf_in, dSdt, composition, collisions, dt, @@ -503,11 +504,10 @@ Function for computing the collision operator ```math \\sum_{s^\\prime} C[F_{s},F_{s^\\prime}] ``` -when -```math -F_{s^\\prime} -``` -is an analytically specified Maxwellian distribution +when \$F_{s^\\prime}\$ +is an analytically specified Maxwellian distribution and +the corresponding Rosenbluth potentials +are specified using analytical results. """ @timeit global_timer fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!( ffs_in, nuref::mk_float, ms::mk_float, Zs::mk_float, @@ -590,11 +590,18 @@ is an analytically specified Maxwellian distribution return nothing end -# solves A x = b for a matrix of the form -# A00 0 A02 -# 0 A11 A12 -# A02 A12 A22 -# appropriate for the moment numerical conserving terms +""" +Function that solves `A x = b` for a matrix of the form +```math +\\begin{array}{ccc} +A_{00} & 0 & A_{02} \\\\ +0 & A_{11} & A_{12} \\\\ +A_{02} & A_{12} & A_{22} \\\\ +\\end{array} +``` +appropriate for the moment numerical conserving terms used in +the Fokker-Planck collision operator. +""" function symmetric_matrix_inverse(A00,A02,A11,A12,A22,b0,b1,b2) # matrix determinant detA = A00*(A11*A22 - A12^2) - A11*A02^2 @@ -615,11 +622,17 @@ function symmetric_matrix_inverse(A00,A02,A11,A12,A22,b0,b1,b2) return x0, x1, x2 end -# solves A x = b for a matrix of the form -# A00 A01 A02 -# A01 A11 A12 -# A02 A12 A22 -# appropriate for the moment numerical conserving terms +""" +Function that solves `A x = b` for a matrix of the form +```math +\\begin{array}{ccc} +A_{00} & A_{01} & A_{02} \\\\ +A_{01} & A_{11} & A_{12} \\\\ +A_{02} & A_{12} & A_{22} \\\\ +\\end{array} +``` +appropriate for moment numerical conserving terms. +""" function symmetric_matrix_inverse(A00,A01,A02,A11,A12,A22,b0,b1,b2) # matrix determinant detA = A00*(A11*A22 - A12^2) - A01*(A01*A22 - A12*A02) + A02*(A01*A12 - A11*A02) @@ -640,6 +653,18 @@ function symmetric_matrix_inverse(A00,A01,A02,A11,A12,A22,b0,b1,b2) return x0, x1, x2 end +""" +Function that applies numerical-error correcting terms to ensure +numerical conservation of the moments `density, upar, pressure` in the self-collision operator. +Modifies the collision operator such that the operator becomes +```math +C_{ss} = C^\\ast_{ss}[F_s,F_{s}] - \\left(x_0 + x_1(v_{\\|}-u_{\\|})+ x_2(v_\\perp^2 +(v_{\\|}-u_{\\|})^2)\\right)F_s +``` +where \$C^\\ast_{ss}[F_s,F_{s}]\$ is the weak-form self-collision operator computed using +the finite-element implementation, \$u_{\\|}\$ is the parallel velocity of \$F_s\$, +and \$x_0,x_1,x_2\$ are parameters that are chosen so that \$C_{ss}\$ +conserves density, parallel velocity and pressure of \$F_s\$. +""" function conserving_corrections!(CC,pdf_in,vpa,vperp,dummy_vpavperp) begin_anyv_region() x0, x1, x2, upar = 0.0, 0.0, 0.0, 0.0 @@ -687,6 +712,15 @@ function conserving_corrections!(CC,pdf_in,vpa,vperp,dummy_vpavperp) end end +""" +Function that applies a numerical-error correcting term to ensure +numerical conservation of the `density` in the collision operator. +```math +C_{ss^\\prime} = C^\\ast_{ss}[F_s,F_{s^\\prime}] - x_0 F_s +``` +where \$C^\\ast_{ss}[F_s,F_{s^\\prime}]\$ is the weak-form collision operator computed using +the finite-element implementation. +""" function density_conserving_correction!(CC,pdf_in,vpa,vperp,dummy_vpavperp) begin_anyv_region() x0 = 0.0 @@ -734,7 +768,7 @@ end """ -allocate the required ancilliary arrays +Function that allocates the required ancilliary arrays for direct integration routines. """ function allocate_fokkerplanck_arrays_direct_integration(vperp,vpa) nvpa = vpa.n @@ -784,7 +818,7 @@ function allocate_fokkerplanck_arrays_direct_integration(vperp,vpa) end """ -function that initialises the arrays needed to calculate the Rosenbluth potentials +Function that initialises the arrays needed to calculate the Rosenbluth potentials by direct integration. As this function is only supported to keep the testing of the direct integration method, the struct 'fka' created here does not contain all of the arrays necessary to compute the weak-form operator. This functionality From c8e7c9b717b05683d5f2d9080a32edc8878f5de4 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 25 Oct 2024 12:59:12 +0100 Subject: [PATCH 2/6] Catch bug where scale and shift factors were calculated with functions rather than using coord.element_scale and coord.element_shift for get_KJ_local!() matrix. Introduce initial explanatory documentation. --- moment_kinetics/src/gauss_legendre.jl | 165 ++++++++++++++++++++------ 1 file changed, 132 insertions(+), 33 deletions(-) diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl index 97c31d54e..34548845a 100644 --- a/moment_kinetics/src/gauss_legendre.jl +++ b/moment_kinetics/src/gauss_legendre.jl @@ -34,9 +34,11 @@ using ..lagrange_polynomials: lagrange_poly_optimised using ..moment_kinetics_structs: weak_discretization_info +#structs for passing around matrices for taking +#the derivatives on Gauss-Legendre points in 1D """ -structs for passing around matrices for taking -the derivatives on Gauss-Legendre points in 1D +A struct for passing around elemental matrices +on Gauss-Legendre points in 1D """ struct gausslegendre_base_info # elementwise differentiation matrix (ngrid*ngrid) @@ -83,6 +85,11 @@ struct gausslegendre_base_info Y31::Array{mk_float,3} end +""" +A struct for Gauss-Legendre arrays needed for global operations in 1D, +contains the struct of elemental matrices for Lobatto and Radau points, +as well as some assembled 1D global matrices. +""" struct gausslegendre_info{TSparse, TSparseCSR, TLU, TLmat, TLmatLU} <: weak_discretization_info lobatto::gausslegendre_base_info radau::gausslegendre_base_info @@ -114,6 +121,9 @@ struct gausslegendre_info{TSparse, TSparseCSR, TLU, TLmat, TLmatLU} <: weak_disc Qmat::Array{mk_float,2} end +""" +Function to create `gausslegendre_info` struct. +""" function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true) lobatto = setup_gausslegendre_pseudospectral_lobatto(coord,collision_operator_dim=collision_operator_dim) radau = setup_gausslegendre_pseudospectral_radau(coord,collision_operator_dim=collision_operator_dim) @@ -153,6 +163,9 @@ function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true) mass_matrix_lu,L_matrix_lu,Qmat) end +""" +Function that fills and `n` x `n` array with the values of the identity matrix `I`. +""" function identity_matrix!(I,n) @. I[:,:] = 0.0 for i in 1:n @@ -161,6 +174,11 @@ function identity_matrix!(I,n) return nothing end +""" +Function that creates the `gausslegendre_base_info` struct for Lobatto points. +If `collision_operator_dim = true`, assign the elemental matrices used to +implement the Fokker-Planck collision operator. +""" function setup_gausslegendre_pseudospectral_lobatto(coord; collision_operator_dim=true) x, w = gausslobatto(coord.ngrid) Dmat = allocate_float(coord.ngrid, coord.ngrid) @@ -231,6 +249,11 @@ function setup_gausslegendre_pseudospectral_lobatto(coord; collision_operator_di K0,K1,K2,P0,P1,P2,D0,Y00,Y01,Y10,Y11,Y20,Y21,Y30,Y31) end +""" +Function that creates the `gausslegendre_base_info` struct for Lobatto points. +If `collision_operator_dim = true`, assign the elemental matrices used to +implement the Fokker-Planck collision operator. +""" function setup_gausslegendre_pseudospectral_radau(coord; collision_operator_dim=true) # Gauss-Radau points on [-1,1) x, w = gaussradau(coord.ngrid) @@ -303,6 +326,10 @@ function setup_gausslegendre_pseudospectral_radau(coord; collision_operator_dim= K0,K1,K2,P0,P1,P2,D0,Y00,Y01,Y10,Y11,Y20,Y21,Y30,Y31) end +""" +A function that takes the first derivative in each element of `coord.grid`, +leaving the result (element-wise) in `coord.scratch_2d`. +""" function elementwise_derivative!(coord, ff, gausslegendre::gausslegendre_info) df = coord.scratch_2d # define local variable nelement for convenience @@ -341,11 +368,19 @@ function elementwise_derivative!(coord, ff, gausslegendre::gausslegendre_info) return nothing end +""" +Wrapper function for element-wise derivatives with advection. +Note that Gauss-Legendre spectral the element method implemented here +does not use upwinding within an element. +""" # Spectral element method does not use upwinding within an element function elementwise_derivative!(coord, ff, adv_fac, spectral::gausslegendre_info) return elementwise_derivative!(coord, ff, spectral) end +""" +Function to perform interpolation on a single element. +""" function single_element_interpolate!(result, newgrid, f, imin, imax, ielement, coord, gausslegendre::gausslegendre_base_info) n_new = length(newgrid) @@ -369,6 +404,9 @@ function single_element_interpolate!(result, newgrid, f, imin, imax, ielement, c return nothing end +""" +Function to carry out a 1D (global) mass matrix solve. +""" function mass_matrix_solve!(f, b, spectral::gausslegendre_info) # invert mass matrix system y = spectral.mass_matrix_lu \ b @@ -377,13 +415,13 @@ function mass_matrix_solve!(f, b, spectral::gausslegendre_info) end """ -Formula for differentiation matrix taken from p196 of Chpt `The Spectral Elemtent Method' of +Formula for Gauss-Legendre-Lobatto differentiation matrix taken from p196 of Chpt `The Spectral Elemtent Method' of `Computational Seismology'. Heiner Igel First Edition. Published in 2017 by Oxford University Press. Or https://doc.nektar.info/tutorials/latest/fundamentals/differentiation/fundamentals-differentiationch2.html -D -- differentiation matrix -x -- Gauss-Legendre-Lobatto points in [-1,1] -ngrid -- number of points per element (incl. boundary points) + D -- differentiation matrix + x -- Gauss-Legendre-Lobatto points in [-1,1] + ngrid -- number of points per element (incl. boundary points) Note that D has does not include a scaling factor """ @@ -408,13 +446,14 @@ function gausslobattolegendre_differentiation_matrix!(D::Array{Float64,2},x::Arr end return nothing end + """ -From +Formula for Gauss-Legendre-Radau differentiation matrix taken from https://doc.nektar.info/tutorials/latest/fundamentals/differentiation/fundamentals-differentiationch2.html -D -- differentiation matrix -x -- Gauss-Legendre-Radau points in [-1,1) -ngrid -- number of points per element (incl. boundary points) + D -- differentiation matrix + x -- Gauss-Legendre-Radau points in [-1,1) + ngrid -- number of points per element (incl. boundary points) Note that D has does not include a scaling factor """ @@ -449,12 +488,12 @@ function gaussradaulegendre_differentiation_matrix!(D::Array{Float64,2},x::Array end """ -Gauss-Legendre derivative at arbitrary x values, for boundary condition on radau points -D0 -- the vector -xj -- the x location where the derivative is evaluated -ngrid -- number of points in x -x -- the grid from -1, 1 -Note that D0 is not scaled to the physical grid +Gauss-Legendre derivative at arbitrary x values, for boundary condition on Radau points. + D0 -- the vector + xj -- the x location where the derivative is evaluated + ngrid -- number of points in x + x -- the grid from -1, 1 +Note that D0 is not scaled to the physical grid with a scaling factor. """ function GaussLegendre_derivative_vector!(D0,xj,ngrid,x,wgts;radau=false) # coefficient in expansion of @@ -482,7 +521,7 @@ function GaussLegendre_derivative_vector!(D0,xj,ngrid,x,wgts;radau=false) end """ -result of the inner product of Legendre polys of order k +Result of the inner product of Legendre polynomials of order k. """ function Legendre_h_n(k) h_n = 2.0/(2.0*k + 1) @@ -491,8 +530,42 @@ end """ -assign abitrary weak inner product matrix Q on a 1D line with Jacobian = 1 -matrix Q acts on a single vector x such that y = Q * x is also a vector +Assign abitrary weak inner product matrix `Q` on a 1D line with Jacobian equal to 1 +matrix `Q` acts on a single vector `x` such that `y = Q * x` is also a vector. + +We use a projection onto Gauss-Legendre polynomials to carry out the calculation +in two steps (see, e.g, S. A. Teukolsky, Short note on the mass matrix for Gauss–Lobatto grid +points, J. Comput. Phys. 283 (2015) 408–413. https://doi.org/10.1016/j.jcp.2014.12.012). +First, we write the desired matrix elements in terms of Legendre polynomials +```math + l_i(x) = \\sum_j \\frac{P_j(x)P_j(x_i)w_i}{\\gamma_j} +``` +with \$w_i\$ the weights from an integration on the Gauss-Legendre-Lobatto (or Radau) points \$x_i\$, +i.e., +```math + \\int^1_{-1} f(x) d x = \\sum_{i} f(x_i)w_i, +``` +and \$\\gamma_j = \\sum_k w_k P_j(x_k)P_j(x_k)\$ the numerical inner-product. +Then, a matrix element can be expressed in integrals over Legendre polynomials +rather than Lagrange polynomials, i.e., +```math + M_{ij} = \\int^1_{-1} l_i(x)l_j(x) d x = \\sum_{mn} \\frac{w_m P_m(x_i) w_n P_n(x_j)}{\\gamma_m\\gamma_n} \\int^{1}_{-1} P_m(x)P_n(x) d x. +``` +Defining +```math + A_{mn} = \\int^{1}_{-1} P_m(x)P_n(x) d x, +``` +we can thus write +```math + M_{ij} = \\sum_{mn} \\frac{w_m P_m(x_i) w_n P_n(x_j)}{\\gamma_m\\gamma_n} A_{mn}. +``` +We can use a quadrature which yields exact results (to machine precision) +to evaluate \$A_{mn}\$ using fast library functions for the Legendre polynomials, +and then carry out the sum \$\\sum_{mn}\$ to obtain exact results (to machine-precision). +Here we use a Gauss-Legendre integration quadrature with exact results up to +polynomials with order \$k_{max} = 4N +1\$, with \$N=\$`ngrid` and the highest order polynomial product +that we integrate is \$P_{N-1}(x)P_{N-1}(x)x^2\$, which has order \$k=2N < k_{max}\$. + """ function GaussLegendre_weak_product_matrix!(QQ::Array{mk_float,2},ngrid,x,wgts,option;radau=false) # coefficient in expansion of @@ -627,9 +700,10 @@ function GaussLegendre_weak_product_matrix!(QQ::Array{mk_float,2},ngrid,x,wgts,o end """ -assign abitrary weak inner product matrix Q on a 1D line with Jacobian = 1 -matrix Q acts on two vectors x1 and x2 such that the quadratic form -y = x1 * Q * x2 is also a vector +Assign abitrary weak (nonlinear) inner product matrix `Q` on a 1D line with Jacobian equal to 1. +matrix `Q` acts on two vectors `x1` and `x2` such that the quadratic form +`y = x1 * Q * x2` is also a vector. See documentation of corresponding function +for linear inner product matrices. """ function GaussLegendre_weak_product_matrix!(QQ::Array{mk_float,3},ngrid,x,wgts,option;radau=false) # coefficient in expansion of @@ -756,10 +830,18 @@ function GaussLegendre_weak_product_matrix!(QQ::Array{mk_float,3},ngrid,x,wgts,o return nothing end +""" +Function for computing the scale factor on a grid with uniformed spaced element boundaries. +Unused. +""" function scale_factor_func(L,nelement_global) return 0.5*L/float(nelement_global) end +""" +Function for computing the shift factor on a grid with uniformed spaced element boundaries. +Unused. +""" function shift_factor_func(L,nelement_global,nelement_local,irank,ielement_local) #ielement_global = ielement_local # for testing + irank*nelement_local ielement_global = ielement_local + irank*nelement_local # proper line for future distributed memory MPI use @@ -767,13 +849,17 @@ function shift_factor_func(L,nelement_global,nelement_local,irank,ielement_local return shift end +""" +Function for finding the elemental index in the global distributed-memory grid. +Distributed-memory for global finite-element operators is not yet supported. +""" function ielement_global_func(nelement_local,irank,ielement_local) return ielement_global = ielement_local + irank*nelement_local end """ -function for setting up the full Gauss-Legendre-Lobatto -grid and collocation point weights +Function for setting up the full Gauss-Legendre-Lobatto +grid and collocation point weights. """ function scaled_gauss_legendre_lobatto_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax) # get Gauss-Legendre-Lobatto points and weights on [-1,1] @@ -802,9 +888,8 @@ function scaled_gauss_legendre_lobatto_grid(ngrid, nelement_local, n_local, elem end """ -function for setting up the full Gauss-Legendre-Radau -grid and collocation point weights -see comments of Gauss-Legendre-Lobatto routine above +Function for setting up the full Gauss-Legendre-Radau +grid and collocation point weights. """ function scaled_gauss_legendre_radau_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax, irank) # get Gauss-Legendre-Lobatto points and weights on [-1,1] @@ -885,7 +970,7 @@ or `dirichlet_bc = true`, `b[1] = f[1]` (except for cylindrical coordinates), `b[end] = f[end]` in the function call, and create new matrices for this purpose -in the gausslegendre_info struct. Currently the Laplacian matrix +in the `gausslegendre_info` struct. Currently the Laplacian matrix is supported with boundary conditions. """ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2}, @@ -964,7 +1049,7 @@ the operators constructed from this function can only be used for differentiation, and not solving 1D ODEs. The shared points in the element assembly are averaged (instead of simply added) to be consistent with the -derivative_elements_to_full_grid!() function in calculus.jl. +`derivative_elements_to_full_grid!()` function in `calculus.jl`. """ function setup_global_strong_form_matrix!(QQ_global::Array{mk_float,2}, lobatto::gausslegendre_base_info, @@ -1024,6 +1109,10 @@ function setup_global_strong_form_matrix!(QQ_global::Array{mk_float,2}, return nothing end +""" +Construction function to provide the appropriate elemental +matrix `Q` to the global matrix assembly functions. +""" function get_QQ_local!(QQ::Array{mk_float,2},ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1057,6 +1146,17 @@ function get_QQ_local!(QQ::Array{mk_float,2},ielement, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `MM` on the \$i^{th}\$ element is +```math + M_{jk} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\varphi_k(v_\\perp) v_\\perp d v_\\perp = \\int^1_{-1} (c_i + x s_i)l_j(x)l_k(x) s_i d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `MM` is +```math + M_{jk} = \\int^{v_\\|^U}_{v_\\|^L} \\varphi_j(v_\\|)\\varphi_k(v_\\|) d v_\\| = \\int^1_{-1} l_j(x)l_k(x) s_i d x. +``` +""" function get_MM_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1159,8 +1259,8 @@ function get_KJ_local!(QQ,ielement, radau::gausslegendre_base_info, coord) - scale_factor = scale_factor_func(coord.L,coord.nelement_global) - shift_factor = shift_factor_func(coord.L,coord.nelement_global,coord.nelement_local,coord.irank,ielement) + 0.5*coord.L + scale_factor = coord.element_scale[ielement] + shift_factor = coord.element_shift[ielement] if coord.name == "vperp" # assume integrals of form int^infty_0 (.) vperp d vperp # extra scale and shift factors required because of vperp^2 in integral if ielement > 1 || coord.irank > 0 # lobatto points @@ -1331,10 +1431,9 @@ function get_PU_local!(QQ,ielement, end """ -construction function for nonlinear diffusion matrices, only +Construction function for nonlinear diffusion matrices, only used in the assembly of the collision operator """ - function get_QQ_local!(QQ::AbstractArray{mk_float,3}, ielement,lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, From 884b1c7eadbd56a8e2aa5b31df52dddd875449bc Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 25 Oct 2024 15:48:51 +0100 Subject: [PATCH 3/6] Complete autodocs for important functions in gauss_legendre.jl --- moment_kinetics/src/gauss_legendre.jl | 134 +++++++++++++++++++++++--- 1 file changed, 123 insertions(+), 11 deletions(-) diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl index 34548845a..c5e5ab677 100644 --- a/moment_kinetics/src/gauss_legendre.jl +++ b/moment_kinetics/src/gauss_legendre.jl @@ -1208,6 +1208,17 @@ function get_SS_local!(QQ,ielement, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `KK` on the \$i^{th}\$ element is +```math + K_{jk} = -\\int^{v_\\perp^U}_{v_\\perp^L} \\left(v_\\perp \\frac{\\partial\\varphi_j(v_\\perp)}{\\partial v_\\perp} + \\varphi_j(v_\\perp) \\right) + \\frac{\\partial\\varphi_k(v_\\perp)}{\\partial v_\\perp} d v_\\perp + = -\\int^1_{-1} ((c_i + x s_i)l_j^\\prime(x) + l_j(x))l_k^\\prime(x) d x /s_i +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `KK` is the same as `LL` (see `get_LL_local!()). +If `explicit_BC_terms = true`, boundary terms arising from integration by parts are included at the extreme boundary points. +""" function get_KK_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1252,8 +1263,15 @@ function get_KK_local!(QQ,ielement, return nothing end -# second derivative matrix with vperp^2 Jacobian factor if -# coord is vperp. Not useful for the vpa coordinate +""" +If called for `coord.name = vperp` elemental matrix `KJ` on the \$i^{th}\$ element is +```math + (KJ)_{jk} = -\\int^{v_\\perp^U}_{v_\\perp^L} \\frac{\\partial\\varphi_j(v_\\perp)}{\\partial v_\\perp}\\frac{\\partial\\varphi_k(v_\\perp)}{\\partial v_\\perp} v_\\perp^2 d v_\\perp + = -\\int^1_{-1} (c_i + x s_i)^2l_j^\\prime(x)l_k^\\prime(x) d x /s_i +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `KJ` is the same as `LL` (see `get_LL_local()!`). +""" function get_KJ_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1278,6 +1296,20 @@ function get_KJ_local!(QQ,ielement, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `LL` on the \$i^{th}\$ element is +```math + L_{jk} = -\\int^{v_\\perp^U}_{v_\\perp^L} \\frac{\\partial\\varphi_j(v_\\perp)}{\\partial v_\\perp}\\frac{\\partial\\varphi_k(v_\\perp)}{\\partial v_\\perp} v_\\perp d v_\\perp + = -\\int^1_{-1} (c_i + x s_i)l_j^\\prime(x)l_k^\\prime(x) d x /s_i +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `LL` is +```math + L_{jk} = -\\int^{v_\\|^U}_{v_\\|^L} \\frac{\\partial\\varphi_j(v_\\|)}{\\partial v_\\|}\\frac{\\partial\\varphi_k(v_\\|)}{\\partial v_\\|} d v_\\| = + -\\int^1_{-1} l_j^\\prime(x)l_k^\\prime(x) d x /s_i. +``` +If `explicit_BC_terms = true`, boundary terms arising from integration by parts are included at the extreme boundary points. +""" function get_LL_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1334,8 +1366,14 @@ function get_DD_local!(QQ, ielement, lobatto::gausslegendre_base_info, return nothing end -# mass matrix without vperp factor (matrix N) -# only useful for the vperp coordinate +""" +If called for `coord.name = vperp` elemental matrix `MN` on the \$i^{th}\$ element is +```math + (MN)_{jk} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\varphi_k(v_\\perp) d v_\\perp = \\int^1_{-1} l_j(x)l_k(x) s_i d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `MN` is the same as `MM` (see `get_MM_local!()`). +""" function get_MN_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1355,8 +1393,14 @@ function get_MN_local!(QQ,ielement, return nothing end -# mass matrix with vperp^2 factor (matrix R) -# only useful for the vperp coordinate +""" +If called for `coord.name = vperp` elemental matrix `MR` on the \$i^{th}\$ element is +```math + (MR)_{jk} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\varphi_k(v_\\perp) v_\\perp^2 d v_\\perp = \\int^1_{-1} (c_i + s_i x)^2 l_j(x)l_k(x) s_i d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `MR` is the same as `MM` (see `get_MM_local!()`). +""" function get_MR_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1381,8 +1425,18 @@ function get_MR_local!(QQ,ielement, return nothing end -# derivative matrix (matrix P, no integration by parts) -# with vperp Jacobian factor if coord is vperp (matrix P) +""" +If called for `coord.name = vperp` elemental matrix `PP` on the \$i^{th}\$ element is +```math + P_{jk} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\frac{\\partial\\varphi_k(v_\\perp)}{\\partial v_\\perp} v_\\perp d v_\\perp + = \\int^1_{-1} (c_i + x s_i)l_j(x)l_k^\\prime(x) d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `PP` is +```math + P_{jk} = \\int^{v_\\|^U}_{v_\\|^L} \\varphi_j(v_\\|)\\frac{\\partial\\varphi_k(v_\\|)}{\\partial v_\\|} d v_\\| = \\int^1_{-1} l_j(x)l_k^\\prime(x) d x. +``` +""" function get_PP_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1403,9 +1457,15 @@ function get_PP_local!(QQ,ielement, return nothing end -# derivative matrix (matrix P, no integration by parts) -# with vperp^2 Jacobian factor if coord is vperp (matrix U) -# not useful for vpa coordinate +""" +If called for `coord.name = vperp` elemental matrix `PP` on the \$i^{th}\$ element is +```math + (PU)_{jk} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\frac{\\partial\\varphi_k(v_\\perp)}{\\partial v_\\perp} v_\\perp^2 d v_\\perp + = \\int^1_{-1} (c_i + x s_i)^2l_j(x)l_k^\\prime(x) d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `PU` is the same as `PP` see `get_PP_local!()`. +""" function get_PU_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1451,6 +1511,19 @@ function get_QQ_local!(QQ::AbstractArray{mk_float,3}, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `YY0` on the \$i^{th}\$ element is +```math + (YY0)_{jkm} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\varphi_k(v_\\perp)\\varphi_m(v_\\perp) v_\\perp d v_\\perp + = \\int^1_{-1} (c_i + x s_i)l_j(x)l_k(x)l_m(x) s_i d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `YY0` is +```math + (YY0)_{jkm} = \\int^{v_\\|^U}_{v_\\|^L} \\varphi_j(v_\\|)\\varphi_k(v_\\|)\\varphi_m(v_\\|) d v_\\| + = \\int^1_{-1} l_j(x)l_k(x)l_m(x) s_i d x. +``` +""" function get_YY0_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1471,6 +1544,19 @@ function get_YY0_local!(QQ,ielement, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `YY1` on the \$i^{th}\$ element is +```math + (YY1)_{jkm} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\varphi_k(v_\\perp)\\frac{\\partial\\varphi_m(v_\\perp)}{\\partial v_\\perp} v_\\perp d v_\\perp + = \\int^1_{-1} (c_i + x s_i)l_j(x)l_k(x)l_m^\\prime(x) d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `YY1` is +```math + (YY1)_{jkm} = \\int^{v_\\|^U}_{v_\\|^L} \\varphi_j(v_\\|)\\varphi_k(v_\\|)\\frac{\\partial\\varphi_m(v_\\|)}{\\partial v_\\|} d v_\\| + = \\int^1_{-1} l_j(x)l_k(x)l_m^\\prime(x) d x. +``` +""" function get_YY1_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1491,6 +1577,19 @@ function get_YY1_local!(QQ,ielement, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `YY2` on the \$i^{th}\$ element is +```math + (YY2)_{jkm} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\frac{\\partial\\varphi_k(v_\\|)}{\\partial v_\\|}\\frac{\\partial\\varphi_m(v_\\perp)}{\\partial v_\\perp} v_\\perp d v_\\perp + = \\int^1_{-1} (c_i + x s_i)l_j(x)l_k^\\prime(x)l_m^\\prime(x) d x/s_i +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `YY2` is +```math + (YY2)_{jkm} = \\int^{v_\\|^U}_{v_\\|^L} \\varphi_j(v_\\|)\\frac{\\partial\\varphi_k(v_\\|)}{\\partial v_\\|}\\frac{\\partial\\varphi_m(v_\\|)}{\\partial v_\\|} d v_\\| + = \\int^1_{-1} l_j(x)l_k^\\prime(x)l_m^\\prime(x) d x /s_i. +``` +""" function get_YY2_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, @@ -1511,6 +1610,19 @@ function get_YY2_local!(QQ,ielement, return nothing end +""" +If called for `coord.name = vperp` elemental matrix `YY3` on the \$i^{th}\$ element is +```math + (YY3)_{jkm} = \\int^{v_\\perp^U}_{v_\\perp^L} \\varphi_j(v_\\perp)\\frac{\\partial\\varphi_k(v_\\|)}{\\partial v_\\|}\\varphi_m(v_\\perp) v_\\perp d v_\\perp + = \\int^1_{-1} (c_i + x s_i)l_j(x)l_k^\\prime(x)l_m(x) d x +``` +with \$c_i\$ and \$s_i\$ the appropriate shift and scale factors, respectively. +Otherwise, if called for any other coordinate elemental matrix `YY3` is +```math + (YY3)_{jkm} = \\int^{v_\\|^U}_{v_\\|^L} \\varphi_j(v_\\|)\\frac{\\partial\\varphi_k(v_\\|)}{\\partial v_\\|}\\varphi_m(v_\\|) d v_\\| + = \\int^1_{-1} l_j(x)l_k^\\prime(x)l_m(x) d x. +``` +""" function get_YY3_local!(QQ,ielement, lobatto::gausslegendre_base_info, radau::gausslegendre_base_info, From 8164762a690ffc804450271a19d1983ecc223c6e Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Sat, 26 Oct 2024 10:59:14 +0100 Subject: [PATCH 4/6] Improve autodocs for fokker_planck_test.jl. --- moment_kinetics/src/fokker_planck_test.jl | 158 ++++++++++++++++++++-- 1 file changed, 145 insertions(+), 13 deletions(-) diff --git a/moment_kinetics/src/fokker_planck_test.jl b/moment_kinetics/src/fokker_planck_test.jl index 103acdc86..9b122a868 100644 --- a/moment_kinetics/src/fokker_planck_test.jl +++ b/moment_kinetics/src/fokker_planck_test.jl @@ -1,7 +1,7 @@ """ -module for including functions used +Module for including functions used in testing the implementation of the -the Full-F Fokker-Planck Collision Operator +the full-F Fokker-Planck collision operator. """ module fokker_planck_test @@ -27,9 +27,18 @@ using ..velocity_moments: get_density # of the Rosenbluth potentials for a shifted Maxwellian # or provide an estimate for collisional coefficients -# G (defined by Del^4 G = -(8/sqrt(pi))*F -# with F = cref^3 pi^(3/2) F_Maxwellian / nref -# the normalised Maxwellian +""" +Function computing G, defined by +```math +\\nabla^4 G = -\\frac{8}{\\sqrt{\\pi}} F +``` +with +```math +F = c_{\\rm ref}^3 \\pi^{3/2} F_{\\rm Maxwellian} / n_{\\rm ref} +``` +the normalised Maxwellian. +See Plasma Confinement, R. D. Hazeltine & J. D. Meiss, 2003, Dover Publications, pg 184, Chpt 5.2, Eqn (5.49). +""" function G_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) # speed variable eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -43,9 +52,18 @@ function G_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) return G*dens*vth end -# H (defined by Del^2 H = -(4/sqrt(pi))*F -# with F = cref^3 pi^(3/2) F_Maxwellian / nref -# the normalised Maxwellian +""" +Function computing H, defined by +```math +\\nabla^2 H = -\\frac{4}{\\sqrt{\\pi}} F +``` +with +```math +F = c_{\\rm ref}^3 \\pi^{3/2} F_{\\rm Maxwellian} / n_{\\rm ref} +``` +the normalised Maxwellian. +See Plasma Confinement, R. D. Hazeltine & J. D. Meiss, 2003, Dover Publications, pg 184, Chpt 5.2, Eqn (5.49). +""" function H_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) # speed variable eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -85,13 +103,27 @@ function dHdeta(eta::mk_float) return dHdeta_fac end -# functions of vpa & vperp +""" +Function computing the normalised speed variable +```math +\\eta = \\frac{\\sqrt{(v_\\| - u_\\|)^2 + v_\\perp^2}}{v_{\\rm th}} +``` +with \$v_{\\rm th} = \\sqrt{2 p / n m}\$ the thermal speed, and \$p\$ the pressure, + \$n\$ the density and \$m\$ the mass. +""" function eta_func(upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) speed = sqrt( (vpa.grid[ivpa] - upar)^2 + vperp.grid[ivperp]^2)/vth return speed end +""" +Function computing +```math +\\frac{\\partial^2 G }{ \\partial v_\\|^2} +``` + for Maxwellian input. See `G_Maxwellian()`. +""" function d2Gdvpa2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -100,6 +132,13 @@ function d2Gdvpa2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return d2Gdvpa2_fac end +""" +Function computing +```math +\\frac{\\partial^2 G}{\\partial v_\\perp \\partial v_\\|} +``` +for Maxwellian input. See `G_Maxwellian()`. +""" function d2Gdvperpdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -108,6 +147,13 @@ function d2Gdvperpdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return d2Gdvperpdvpa_fac end +""" +Function computing +```math +\\frac{\\partial^2 G}{\\partial v_\\perp^2} +``` +for Maxwellian input. See `G_Maxwellian()`. +""" function d2Gdvperp2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -116,6 +162,13 @@ function d2Gdvperp2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return d2Gdvperp2_fac end +""" +Function computing +```math +\\frac{\\partial G}{\\partial v_\\perp} +``` +for Maxwellian input. See `G_Maxwellian()`. +""" function dGdvperp_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -123,6 +176,13 @@ function dGdvperp_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial H}{\\partial v_\\perp} +``` +for Maxwellian input. See `H_Maxwellian()`. +""" function dHdvperp_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -130,6 +190,13 @@ function dHdvperp_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial H}{\\partial v_\\|} +``` +for Maxwellian input. See `H_Maxwellian()`. +""" function dHdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -137,6 +204,9 @@ function dHdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing \$ F_{\\rm Maxwellian} \$. +""" function F_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -144,6 +214,13 @@ function F_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial F}{\\partial v_\\|} +``` +for \$ F = F_{\\rm Maxwellian}\$. +""" function dFdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -151,6 +228,13 @@ function dFdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial F}{\\partial v_\\perp} +``` +for \$ F = F_{\\rm Maxwellian}\$. +""" function dFdvperp_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -158,6 +242,13 @@ function dFdvperp_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial^2 F}{\\partial v_\\perp \\partial v_\\|} +``` +for \$ F = F_{\\rm Maxwellian}\$. +""" function d2Fdvperpdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -165,6 +256,13 @@ function d2Fdvperpdvpa_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial^2 F}{\\partial v_\\|^2} +``` +for \$ F = F_{\\rm Maxwellian}\$. +""" function d2Fdvpa2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -172,6 +270,13 @@ function d2Fdvpa2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Function computing +```math +\\frac{\\partial^2 F}{\\partial v_\\perp^2}. +``` +for \$ F = F_{\\rm Maxwellian}\$. +""" function d2Fdvperp2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, vpa,vperp,ivpa,ivperp) eta = eta_func(upar,vth,vpa,vperp,ivpa,ivperp) @@ -179,6 +284,10 @@ function d2Fdvperp2_Maxwellian(dens::mk_float,upar::mk_float,vth::mk_float, return fac end +""" +Calculates the fully expanded form of the collision operator \$C_{s s^\\prime}[F_s,F_{s^\\prime}]\$ given Maxwellian input \$F_s\$ and \$F_{s^\\prime}\$. +The input Maxwellians are specified through their moments. +""" function Cssp_Maxwellian_inputs(denss::mk_float,upars::mk_float,vths::mk_float,ms::mk_float, denssp::mk_float,uparsp::mk_float,vthsp::mk_float,msp::mk_float, nussp::mk_float,vpa,vperp,ivpa,ivperp) @@ -210,6 +319,10 @@ function Cssp_Maxwellian_inputs(denss::mk_float,upars::mk_float,vths::mk_float,m return Cssp_Maxwellian end +""" +Calculates the collisional flux \$\\Gamma_\\|\$ given Maxwellian input \$F_s\$ and \$F_{s^\\prime}\$. +The input Maxwellians are specified through their moments. +""" function Cflux_vpa_Maxwellian_inputs(ms::mk_float,denss::mk_float,upars::mk_float,vths::mk_float, msp::mk_float,denssp::mk_float,uparsp::mk_float,vthsp::mk_float, vpa,vperp,ivpa,ivperp) @@ -224,6 +337,10 @@ function Cflux_vpa_Maxwellian_inputs(ms::mk_float,denss::mk_float,upars::mk_floa return Cflux end +""" +Calculates the collisional flux \$\\Gamma_\\perp\$ given Maxwellian input \$F_s\$ and \$F_{s^\\prime}\$. +The input Maxwellians are specified through their moments. +""" function Cflux_vperp_Maxwellian_inputs(ms::mk_float,denss::mk_float,upars::mk_float,vths::mk_float, msp::mk_float,denssp::mk_float,uparsp::mk_float,vthsp::mk_float, vpa,vperp,ivpa,ivperp) @@ -240,7 +357,8 @@ end """ Function calculating the fully expanded form of the collision operator -taking floats as arguments. This function is designed to be used at the +taking as arguments the derivatives of \$F_s\$, \$G_{s^\\prime}\$ and \$H_{s^\\prime}\$. +This function is designed to be used at the lowest level of a coordinate loop, with derivatives and integrals all previously calculated. """ @@ -259,7 +377,7 @@ end """ -calculates the collisional fluxes given input F_s and G_sp, H_sp +Calculates the collisional fluxes given input \$F_s\$ and \$G_{s^\\prime}\$, \$H_{s^\\prime}\$. """ function calculate_collisional_fluxes(F,dFdvpa,dFdvperp, d2Gdvpa2,d2Gdvperpdvpa,d2Gdvperp2,dHdvpa,dHdvperp, @@ -273,10 +391,11 @@ function calculate_collisional_fluxes(F,dFdvpa,dFdvperp, end +# Below are functions which are used for storing and printing data from the tests + """ -Below are functions which are used for storing and printing data from the tests +Function to print the maximum error \${\\rm MAX}(|f_{\\rm numerical}-f_{\\rm exact}|)\$. """ - function print_test_data(func_exact,func_num,func_err,func_name) @. func_err = abs(func_num - func_exact) max_err = maximum(func_err) @@ -284,6 +403,13 @@ function print_test_data(func_exact,func_num,func_err,func_name) return max_err end +""" +Function to print the maximum error \${\\rm MAX}(|f_{\\rm numerical}-f_{\\rm exact}|)\$ and the +\$L_2\$ norm of the error +```math +\\sqrt{\\int (f - f_{\\rm exact})^2 v_\\perp d v_\\perp d v_\\|/\\int v_\\perp d v_\\perp d v_\\|}. +``` +""" function print_test_data(func_exact,func_num,func_err,func_name,vpa,vperp,dummy;print_to_screen=true) @. func_err = abs(func_num - func_exact) max_err = maximum(func_err) @@ -340,6 +466,9 @@ function allocate_error_data() moments) end +""" +Utility function that saves error data to a HDF5 file for later use. +""" function save_fkpl_error_data(outdir,ncore,ngrid,nelement_list, max_C_err, max_H_err, max_G_err, max_dHdvpa_err, max_dHdvperp_err, max_d2Gdvperp2_err, max_d2Gdvpa2_err, max_d2Gdvperpdvpa_err, max_dGdvperp_err, @@ -384,6 +513,9 @@ function save_fkpl_error_data(outdir,ncore,ngrid,nelement_list, return nothing end +""" +Utility function that saves error data to a HDF5 file for later use. +""" function save_fkpl_integration_error_data(outdir,ncore,ngrid,nelement_list, max_dfsdvpa_err, max_dfsdvperp_err, max_d2fsdvperpdvpa_err, max_H_err, max_G_err, max_dHdvpa_err, max_dHdvperp_err, From 0538cf08332fde236f9c273a249ac9e4d1e30bc5 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Sun, 27 Oct 2024 09:04:04 +0000 Subject: [PATCH 5/6] Improve autodocs for fokker_planck_calculus, remove export statement for unused function. --- moment_kinetics/src/fokker_planck.jl | 2 +- moment_kinetics/src/fokker_planck_calculus.jl | 372 ++++++++++++++---- test_scripts/2D_FEM_assembly_test.jl | 2 +- 3 files changed, 291 insertions(+), 85 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 521dca83d..9a19affb9 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -66,7 +66,7 @@ using ..fokker_planck_calculus: assemble_explicit_collision_operator_rhs_paralle using ..fokker_planck_calculus: assemble_explicit_collision_operator_rhs_parallel_analytical_inputs! using ..fokker_planck_calculus: calculate_YY_arrays, enforce_vpavperp_BCs! using ..fokker_planck_calculus: calculate_rosenbluth_potential_boundary_data! -using ..fokker_planck_calculus: enforce_zero_bc!, elliptic_solve!, algebraic_solve! +using ..fokker_planck_calculus: elliptic_solve!, algebraic_solve! using ..fokker_planck_calculus: calculate_rosenbluth_potentials_via_elliptic_solve! using ..fokker_planck_test: Cssp_fully_expanded_form, calculate_collisional_fluxes, H_Maxwellian, dGdvperp_Maxwellian using ..fokker_planck_test: d2Gdvpa2_Maxwellian, d2Gdvperpdvpa_Maxwellian, d2Gdvperp2_Maxwellian, dHdvpa_Maxwellian, dHdvperp_Maxwellian diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index f11dcf43c..85939a5fa 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -1,9 +1,9 @@ """ -module for functions used +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`](@ref). +the full-F Fokker-Planck collision operator [`moment_kinetics.fokker_planck`](@ref). Parallelisation of the collision operator uses a special 'anyv' region type, see [Collision operator and `anyv` region](@ref). @@ -25,7 +25,6 @@ export calculate_rosenbluth_potentials_via_elliptic_solve! # testing export calculate_rosenbluth_potential_boundary_data_exact! -export enforce_zero_bc! export allocate_rosenbluth_potential_boundary_data export calculate_rosenbluth_potential_boundary_data_exact! export test_rosenbluth_potential_boundary_data @@ -71,10 +70,11 @@ function print_vector(vector,name::String,m::mk_int) end """ -a struct of dummy arrays and precalculated coefficients -for the strong-form Fokker-Planck collision operator +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 fokkerplanck_arrays_direct_integration_struct G0_weights::MPISharedArray{mk_float,4} G1_weights::MPISharedArray{mk_float,4} @@ -103,8 +103,8 @@ struct fokkerplanck_arrays_direct_integration_struct end """ -a struct to contain the integration weights for the boundary points -in the (vpa,vperp) domain +Struct to contain the integration weights for the boundary points +in the `(vpa,vperp)` domain. """ struct boundary_integration_weights_struct lower_vpa_boundary::MPISharedArray{mk_float,3} @@ -113,8 +113,8 @@ struct boundary_integration_weights_struct end """ -a struct used for calculating the integration weights for the -boundary of the velocity space domain in (vpa,vperp) coordinates +Struct used for storing the integration weights for the +boundary of the velocity space domain in `(vpa,vperp)` coordinates. """ struct fokkerplanck_boundary_data_arrays_struct G0_weights::boundary_integration_weights_struct @@ -128,12 +128,20 @@ struct fokkerplanck_boundary_data_arrays_struct dfdvperp::MPISharedArray{mk_float,2} end +""" +Struct to store the `(vpa,vperp)` boundary data for an +individual Rosenbluth potential. +""" struct vpa_vperp_boundary_data lower_boundary_vpa::MPISharedArray{mk_float,1} upper_boundary_vpa::MPISharedArray{mk_float,1} upper_boundary_vperp::MPISharedArray{mk_float,1} end +""" +Struct to store the boundary data for all of the +Rosenbluth potentials required for the calculation. +""" struct rosenbluth_potential_boundary_data H_data::vpa_vperp_boundary_data dHdvpa_data::vpa_vperp_boundary_data @@ -145,6 +153,15 @@ struct rosenbluth_potential_boundary_data d2Gdvpa2_data::vpa_vperp_boundary_data end +""" +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 YY_collision_operator_arrays # let phi_j(vperp) be the jth Lagrange basis function, # and phi'_j(vperp) the first derivative of the Lagrange basis function @@ -168,8 +185,8 @@ struct YY_collision_operator_arrays end """ -a struct of dummy arrays and precalculated coefficients -for the weak-form Fokker-Planck collision operator +Struct of dummy arrays and precalculated coefficients +for the finite-element weak-form Fokker-Planck collision operator. """ struct fokkerplanck_weakform_arrays_struct{M <: AbstractSparseArray{mk_float,mk_int,N} where N} # boundary weights (Green's function) data @@ -220,6 +237,9 @@ struct fokkerplanck_weakform_arrays_struct{M <: AbstractSparseArray{mk_float,mk_ dFdvperp::MPISharedArray{mk_float,2} end +""" +Function to allocate a `boundary_integration_weights_struct`. +""" function allocate_boundary_integration_weight(vpa,vperp) nvpa = vpa.n nvperp = vperp.n @@ -230,6 +250,9 @@ function allocate_boundary_integration_weight(vpa,vperp) upper_vpa_boundary, upper_vperp_boundary) end +""" +Function to allocate at `fokkerplanck_boundary_data_arrays_struct`. +""" function allocate_boundary_integration_weights(vpa,vperp) G0_weights = allocate_boundary_integration_weight(vpa,vperp) G1_weights = allocate_boundary_integration_weight(vpa,vperp) @@ -260,7 +283,8 @@ end """ -function that precomputes the required integration weights +Function that precomputes the required integration weights in the whole of +`(vpa,vperp)` for the direct integration method of computing the Rosenbluth potentials. """ function init_Rosenbluth_potential_integration_weights!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights,vperp,vpa;print_to_screen=true) @@ -314,9 +338,9 @@ function init_Rosenbluth_potential_integration_weights!(G0_weights,G1_weights,H0 end """ -function for getting the basic quadratures used for the +Function for getting the basic quadratures used for the numerical integration of the Lagrange polynomials and the -Green's function. +integration kernals. """ function setup_basic_quadratures(vpa,vperp;print_to_screen=true) @serial_region begin @@ -341,8 +365,7 @@ end """ -function for getting the indices used to choose the integration -quadrature +Function for getting the indices used to choose the integration quadrature. """ function get_element_limit_indices(ivpa,ivperp,vpa,vperp) nelement_vpa, ngrid_vpa = vpa.nelement_local, vpa.ngrid @@ -359,9 +382,11 @@ function get_element_limit_indices(ivpa,ivperp,vpa,vperp) return igrid_vpa, ielement_vpa, ielement_vpa_low, ielement_vpa_hi, igrid_vperp, ielement_vperp, ielement_vperp_low, ielement_vperp_hi end + """ -function that precomputes the required integration weights -only along the velocity space boundaries +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 init_Rosenbluth_potential_boundary_integration_weights!(G0_weights, G1_weights,H0_weights,H1_weights,H2_weights,H3_weights,vpa,vperp;print_to_screen=true) @@ -508,12 +533,14 @@ function get_nodes(coord,iel) return nodes end -# 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. +""" +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. +""" function get_scaled_x_w_with_divergences!(x_scaled, w_scaled, x_legendre, w_legendre, x_laguerre, w_laguerre, node_min, node_max, nodes, igrid_coord, coord_val) #println("nodes ",nodes) zero = 1.0e-10 @@ -604,9 +631,12 @@ function get_scaled_x_w_with_divergences!(x_scaled, w_scaled, x_legendre, w_lege #println("w_scaled",w_scaled) return nquad_coord end -# 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 grid and integration weights assuming +no divergences of the function on the 1D element. Gauss-Legendre +quadrature is used for the entire element. +""" function get_scaled_x_w_no_divergences!(x_scaled, w_scaled, x_legendre, w_legendre, node_min, node_max) @. x_scaled = 0.0 @. w_scaled = 0.0 @@ -621,30 +651,46 @@ function get_scaled_x_w_no_divergences!(x_scaled, w_scaled, x_legendre, w_legend return nquad end -# function returns 1 if igrid = 1 or 0 if 1 < igrid <= ngrid +""" +Function returns `1` if `igrid = 1` or `0` if `1 < igrid <= ngrid`. +""" function ng_low(igrid,ngrid) return floor(mk_int, (ngrid - igrid)/(ngrid - 1)) end -# function returns 1 if igrid = ngrid or 0 if 1 =< igrid < ngrid + +""" +Function returns `1` if `igrid = ngrid` or `0` if `1 =< igrid < ngrid`. +""" function ng_hi(igrid,ngrid) return floor(mk_int, igrid/ngrid) end -# function returns 1 for nelement >= ielement > 1, 0 for ielement =1 + +""" +Function returns `1` for `nelement >= ielement > 1`, `0` for `ielement = 1`. +""" function nel_low(ielement,nelement) return floor(mk_int, (ielement - 2 + nelement)/nelement) end -# function returns 1 for nelement > ielement >= 1, 0 for ielement =nelement + +""" +Function returns `1` for `nelement > ielement >= 1`, `0` for `ielement = nelement`. +""" function nel_hi(ielement,nelement) return 1- floor(mk_int, ielement/nelement) end -# base level function for computing the Green's function weights -# note the definitions of ellipe & ellipk -# `https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.ellipe` -# `https://specialfunctions.juliamath.org/stable/functions_list/#SpecialFunctions.ellipk` -# `ellipe(m) = \int^{\pi/2}\_0 \sqrt{ 1 - m \sin^2(\theta)} d \theta` -# `ellipe(k) = \int^{\pi/2}\_0 \frac{1}{\sqrt{ 1 - m \sin^2(\theta)}} d \theta` - +""" +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` +```math +E(m) = \\int^{\\pi/2}_0 \\sqrt{ 1 - m \\sin^2(\\theta)} d \\theta +``` +```math +K(m) = \\int^{\\pi/2}_0 \\frac{1}{\\sqrt{ 1 - m \\sin^2(\\theta)}} d \\theta +``` +""" function local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, nquad_vpa,ielement_vpa,vpa, # info about primed vpa grids nquad_vperp,ielement_vperp,vperp, # info about primed vperp grids @@ -733,6 +779,13 @@ function local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights, return nothing end +""" +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 loop_over_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, vpa,ielement_vpa_low,ielement_vpa_hi, # info about primed vperp grids vperp,ielement_vperpp, # info about primed vperp grids @@ -784,6 +837,11 @@ function loop_over_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weights,H2_ return nothing end +""" +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 loop_over_vpa_elements_no_divergences!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, vpa,ielement_vpa_low,ielement_vpa_hi, # info about primed vperp grids nquad_vperp,ielement_vperpp,vperp_nodes,vperp, # info about primed vperp grids @@ -805,6 +863,15 @@ function loop_over_vpa_elements_no_divergences!(G0_weights,G1_weights,H0_weights return nothing end +""" +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. +""" function loop_over_vperp_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, vpa,ielement_vpa_low,ielement_vpa_hi, # info about primed vpa grids vperp,ielement_vperp_low,ielement_vperp_hi, # info about primed vperp grids @@ -854,12 +921,14 @@ function loop_over_vperp_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weigh return nothing end -# The function loop_over_vperp_vpa_elements_no_divergences!() was 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 the problems with the pure integration method of computing the -# Rosenbluth potentials are understood. +""" +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 loop_over_vperp_vpa_elements_no_divergences!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, vpa,ielement_vpa_low,ielement_vpa_hi, # info about primed vpa grids vperp,ielement_vperp_low,ielement_vperp_hi, # info about primed vperp grids @@ -926,9 +995,11 @@ function ivpa_func(ic::mk_int,nvpa::mk_int) return ivpa end -# function that returns the sparse matrix index -# used to directly construct the nonzero entries -# of a 2D assembled sparse matrix +""" +Function that returns the sparse matrix index +used to directly construct the nonzero entries +of a 2D assembled sparse matrix. +""" function icsc_func(ivpa_local::mk_int,ivpap_local::mk_int, ielement_vpa::mk_int, ngrid_vpa::mk_int,nelement_vpa::mk_int, @@ -946,6 +1017,9 @@ function icsc_func(ivpa_local::mk_int,ivpap_local::mk_int, return icsc end +""" +Struct to contain data needed to create a sparse matrix. +""" struct sparse_matrix_constructor # the Ith row II::Array{mk_float,1} @@ -955,6 +1029,9 @@ struct sparse_matrix_constructor SS::Array{mk_float,1} end +""" +Function to allocate an instance of `sparse_matrix_constructor`. +""" function allocate_sparse_matrix_constructor(nsparse::mk_int) II = Array{mk_int,1}(undef,nsparse) @. II = 0 @@ -965,12 +1042,20 @@ function allocate_sparse_matrix_constructor(nsparse::mk_int) return sparse_matrix_constructor(II,JJ,SS) end +""" +Function to assign data to an instance of `sparse_matrix_constructor`. +""" function assign_constructor_data!(data::sparse_matrix_constructor,icsc::mk_int,ii::mk_int,jj::mk_int,ss::mk_float) data.II[icsc] = ii data.JJ[icsc] = jj data.SS[icsc] = ss return nothing end + +""" +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 assemble_constructor_data!(data::sparse_matrix_constructor,icsc::mk_int,ii::mk_int,jj::mk_int,ss::mk_float) data.II[icsc] = ii data.JJ[icsc] = jj @@ -978,10 +1063,17 @@ function assemble_constructor_data!(data::sparse_matrix_constructor,icsc::mk_int return nothing end +""" +Wrapper function to create a sparse matrix with an instance of `sparse_matrix_constructor` +and `sparse()`. +""" function create_sparse_matrix(data::sparse_matrix_constructor) return sparse(data.II,data.JJ,data.SS) end +""" +Function to allocate an instance of `vpa_vperp_boundary_data`. +""" function allocate_boundary_data(vpa,vperp) # The following velocity-space-sized buffer arrays are used to evaluate the # collision operator for a single species at a single spatial point. They are @@ -997,7 +1089,10 @@ function allocate_boundary_data(vpa,vperp) upper_boundary_vpa,upper_boundary_vperp) end - +""" +Function to assign precomputed (exact) data to an instance +of `vpa_vperp_boundary_data`. Used in testing. +""" function assign_exact_boundary_data!(func_data::vpa_vperp_boundary_data, func_exact,vpa,vperp) begin_anyv_region() @@ -1014,7 +1109,10 @@ function assign_exact_boundary_data!(func_data::vpa_vperp_boundary_data, end return nothing end - + +""" +Function to allocate an instance of `rosenbluth_potential_boundary_data`. +""" function allocate_rosenbluth_potential_boundary_data(vpa,vperp) H_data = allocate_boundary_data(vpa,vperp) dHdvpa_data = allocate_boundary_data(vpa,vperp) @@ -1029,6 +1127,10 @@ function allocate_rosenbluth_potential_boundary_data(vpa,vperp) d2Gdvperpdvpa_data,d2Gdvpa2_data) end +""" +Function to assign data to an instance of `rosenbluth_potential_boundary_data`, in place, +without allocation. Used in testing. +""" function calculate_rosenbluth_potential_boundary_data_exact!(rpbd::rosenbluth_potential_boundary_data, H_exact,dHdvpa_exact,dHdvperp_exact,G_exact,dGdvperp_exact, d2Gdvperp2_exact,d2Gdvperpdvpa_exact,d2Gdvpa2_exact, @@ -1044,7 +1146,13 @@ function calculate_rosenbluth_potential_boundary_data_exact!(rpbd::rosenbluth_po return nothing end - +""" +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 calculate_boundary_data!(func_data::vpa_vperp_boundary_data, weight::MPISharedArray{mk_float,4},func_input,vpa,vperp) nvpa = vpa.n @@ -1073,6 +1181,12 @@ function calculate_boundary_data!(func_data::vpa_vperp_boundary_data, return nothing end +""" +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 calculate_boundary_data!(func_data::vpa_vperp_boundary_data, weight::boundary_integration_weights_struct, func_input,vpa,vperp) @@ -1103,6 +1217,11 @@ function calculate_boundary_data!(func_data::vpa_vperp_boundary_data, return nothing end +""" +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 calculate_rosenbluth_potential_boundary_data!(rpbd::rosenbluth_potential_boundary_data, fkpl::Union{fokkerplanck_arrays_direct_integration_struct,fokkerplanck_boundary_data_arrays_struct},pdf,vpa,vperp,vpa_spectral,vperp_spectral; calculate_GG=false,calculate_dGdvperp=false) @@ -1140,6 +1259,11 @@ function calculate_rosenbluth_potential_boundary_data!(rpbd::rosenbluth_potentia return nothing end +""" +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()`. +""" function test_rosenbluth_potential_boundary_data(rpbd::rosenbluth_potential_boundary_data, rpbd_exact::rosenbluth_potential_boundary_data,vpa,vperp;print_to_screen=true) @@ -1158,6 +1282,10 @@ function test_rosenbluth_potential_boundary_data(rpbd::rosenbluth_potential_boun return max_H_err, max_dHdvpa_err, max_dHdvperp_err, max_G_err, max_dGdvperp_err, max_d2Gdvperp2_err, max_d2Gdvperpdvpa_err, max_d2Gdvpa2_err end +""" +Function to compute the maximum error \${\\rm MAX}|f_{\\rm numerical}-f_{\\rm exact}|\$ for +instances of `vpa_vperp_boundary_data`. +""" function test_boundary_data(func,func_exact,func_name,vpa,vperp,buffer_vpa,buffer_vperp_1,buffer_vperp_2,print_to_screen) nvpa = vpa.n nvperp = vperp.n @@ -1197,6 +1325,10 @@ function get_global_compound_index(vpa,vperp,ielement_vpa,ielement_vperp,ivpa_lo return ic_global end +""" +Unused function. Sets `f(vpa,vperp)` to zero at the boundaries +in `(vpa,vperp)`. +""" function enforce_zero_bc!(fvpavperp,vpa,vperp;impose_BC_at_zero_vperp=false) # lower vpa boundary @loop_vperp ivperp begin @@ -1221,6 +1353,11 @@ function enforce_zero_bc!(fvpavperp,vpa,vperp;impose_BC_at_zero_vperp=false) end end +""" +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. +""" function enforce_dirichlet_bc!(fvpavperp,vpa,vperp,f_bc;dirichlet_vperp_lower_boundary=false) # lower vpa boundary for ivperp ∈ 1:vperp.n @@ -1245,6 +1382,10 @@ function enforce_dirichlet_bc!(fvpavperp,vpa,vperp,f_bc;dirichlet_vperp_lower_bo end end +""" +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 enforce_dirichlet_bc!(fvpavperp,vpa,vperp,f_bc::vpa_vperp_boundary_data) # lower vpa boundary for ivperp ∈ 1:vperp.n @@ -1263,6 +1404,16 @@ function enforce_dirichlet_bc!(fvpavperp,vpa,vperp,f_bc::vpa_vperp_boundary_data return nothing end +""" +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 assemble_matrix_operators_dirichlet_bc(vpa,vperp,vpa_spectral,vperp_spectral;print_to_screen=true) nc_global = vpa.n*vperp.n # Assemble a 2D mass matrix in the global compound coordinate @@ -1486,6 +1637,16 @@ function assemble_matrix_operators_dirichlet_bc(vpa,vperp,vpa_spectral,vperp_spe PPpar2D_sparse, MMparMNperp2D_sparse end +""" +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 assemble_matrix_operators_dirichlet_bc_sparse(vpa,vperp,vpa_spectral,vperp_spectral;print_to_screen=true) # Assemble a 2D mass matrix in the global compound coordinate nc_global = vpa.n*vperp.n @@ -1725,6 +1886,11 @@ function assemble_matrix_operators_dirichlet_bc_sparse(vpa,vperp,vpa_spectral,vp PPpar2D_sparse, MMparMNperp2D_sparse end +""" +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 calculate_YY_arrays(vpa,vperp,vpa_spectral,vperp_spectral) YY0perp = Array{mk_float,4}(undef,vperp.ngrid,vperp.ngrid,vperp.ngrid,vperp.nelement_local) YY1perp = Array{mk_float,4}(undef,vperp.ngrid,vperp.ngrid,vperp.ngrid,vperp.nelement_local) @@ -1752,6 +1918,12 @@ function calculate_YY_arrays(vpa,vperp,vpa_spectral,vperp_spectral) YY0par,YY1par,YY2par,YY3par) end +""" +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 assemble_explicit_collision_operator_rhs_serial!(rhsvpavperp,pdfs,d2Gspdvpa2,d2Gspdvperpdvpa, d2Gspdvperp2,dHspdvpa,dHspdvperp,ms,msp,nussp, vpa,vperp,YY_arrays::YY_collision_operator_arrays) @@ -1812,6 +1984,14 @@ function assemble_explicit_collision_operator_rhs_serial!(rhsvpavperp,pdfs,d2Gsp return nothing end +""" +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 assemble_explicit_collision_operator_rhs_parallel!(rhsvpavperp,pdfs,d2Gspdvpa2,d2Gspdvperpdvpa, d2Gspdvperp2,dHspdvpa,dHspdvperp,ms,msp,nussp, vpa,vperp,YY_arrays::YY_collision_operator_arrays) @@ -1859,6 +2039,9 @@ function assemble_explicit_collision_operator_rhs_parallel!(rhsvpavperp,pdfs,d2G return nothing end +""" +The inner-most loop of the parallel collision operator assembly. Called in `assemble_explicit_collision_operator_rhs_parallel!()`. +""" function assemble_explicit_collision_operator_rhs_parallel_inner_loop( nussp, ms, msp, YY0perp, YY0par, YY1perp, YY1par, YY2perp, YY2par, YY3perp, YY3par, pdfs, d2Gspdvpa2, d2Gspdvperpdvpa, d2Gspdvperp2, dHspdvpa, dHspdvperp, @@ -1898,6 +2081,12 @@ function assemble_explicit_collision_operator_rhs_parallel_inner_loop( return result end +""" +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()`. +""" function assemble_explicit_collision_operator_rhs_parallel_analytical_inputs!(rhsvpavperp,pdfs,dpdfsdvpa,dpdfsdvperp,d2Gspdvpa2,d2Gspdvperpdvpa, d2Gspdvperp2,dHspdvpa,dHspdvperp,ms,msp,nussp, vpa,vperp,YY_arrays::YY_collision_operator_arrays) @@ -1945,6 +2134,9 @@ function assemble_explicit_collision_operator_rhs_parallel_analytical_inputs!(rh return nothing end +""" +The inner-most loop of `assemble_explicit_collision_operator_rhs_parallel_analytical_inputs!()`. +""" # Separate function for inner loop, possible optimization?? function assemble_explicit_collision_operator_rhs_parallel_analytical_inputs_inner_loop( nussp, ms, msp, pdfs, dpdfsdvpa, dpdfsdvperp, d2Gspdvperp2, @@ -1987,18 +2179,20 @@ function assemble_explicit_collision_operator_rhs_parallel_analytical_inputs_inn return result end - -# 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. +""" +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. +""" function elliptic_solve!(field,source,boundary_data::vpa_vperp_boundary_data, lu_object_lhs,matrix_rhs,rhsvpavperp,vpa,vperp) # assemble the rhs of the weak system @@ -2042,14 +2236,16 @@ function elliptic_solve!(field,source_1,source_2,boundary_data::vpa_vperp_bounda return nothing end -# 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. +""" +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 algebraic_solve!(field,source_1,source_2,boundary_data::vpa_vperp_boundary_data, lu_object_lhs,matrix_rhs_1,matrix_rhs_2,rhs,vpa,vperp) @@ -2073,6 +2269,15 @@ function algebraic_solve!(field,source_1,source_2,boundary_data::vpa_vperp_bound return nothing end +""" +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. +""" function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvperp, d2Gdvpa2,dGdvperp,d2Gdvperpdvpa,d2Gdvperp2,ffsp_in, vpa,vperp,vpa_spectral,vperp_spectral,fkpl_arrays::fokkerplanck_weakform_arrays_struct; @@ -2199,7 +2404,8 @@ function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvpe end """ -function to calculate Rosenbluth potentials by direct integration +Function to calculate Rosenbluth potentials in the entire +domain of `(vpa,vperp)` by direct integration. """ function calculate_rosenbluth_potentials_via_direct_integration!(GG,HH,dHdvpa,dHdvperp, @@ -2277,8 +2483,9 @@ function calculate_rosenbluth_integrals!(GG,d2Gspdvpa2,dGspdvperp,d2Gspdvperpdvp end """ -function to enforce boundary conditions on the collision operator -result to be consistent with the boundary conditions imposed on the the pdf +Function to enforce boundary conditions on the collision operator +result to be consistent with the boundary conditions imposed on the +distribution function. """ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral) nvpa = vpa.n @@ -2314,20 +2521,19 @@ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral) end """ -function to interpolate f(vpa,vperp) from one +Function to interpolate `f(vpa,vperp)` from one velocity grid to another, assuming that both -grids are represented by vpa, vperp in normalised units, +grids are represented by `(vpa,vperp)` in normalised units, but have different normalisation factors -defining the meaning of these grids in physical units. +defining the meaning of these grids in physical units. E.g., -E.g. vpai, vperpi = ci * vpa, ci * vperp + 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 +with `ci = sqrt(Ti/mi)`, `ce = sqrt(Te/mi)` +`scalefac = ci / ce` is the ratio of the +two reference speeds. """ function interpolate_2D_vspace!(pdf_out,pdf_in,vpa,vperp,scalefac) @@ -2403,7 +2609,7 @@ end #end """ -function to find the element in which x sits +Function to find the element in which x sits. """ function ielement_loopup(x,coord) xebs = coord.element_boundaries diff --git a/test_scripts/2D_FEM_assembly_test.jl b/test_scripts/2D_FEM_assembly_test.jl index bdb499490..34a19c3cd 100644 --- a/test_scripts/2D_FEM_assembly_test.jl +++ b/test_scripts/2D_FEM_assembly_test.jl @@ -30,7 +30,7 @@ using moment_kinetics.fokker_planck_test: print_test_data, fkpl_error_data, allo using moment_kinetics.fokker_planck_test: save_fkpl_error_data using moment_kinetics.fokker_planck_calculus: elliptic_solve! -using moment_kinetics.fokker_planck_calculus: enforce_zero_bc!, allocate_rosenbluth_potential_boundary_data +using moment_kinetics.fokker_planck_calculus: allocate_rosenbluth_potential_boundary_data using moment_kinetics.fokker_planck_calculus: calculate_rosenbluth_potential_boundary_data!, calculate_rosenbluth_potential_boundary_data_exact! using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary_data, enforce_vpavperp_BCs! using moment_kinetics.fokker_planck_calculus: calculate_rosenbluth_potentials_via_elliptic_solve! From 3801743d7e6883c0dd020d35605e29a55f4db378 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Sun, 27 Oct 2024 09:39:01 +0000 Subject: [PATCH 6/6] Add notes on the Fokker-Planck collision operator, linking to the ExCALIBUR report and giving input parameter information. Modify the index to feature this note and the note on magnetic geometry. --- docs/src/fokker_planck_notes.md | 37 +++++++++++++++++++++++++++++++++ docs/src/index.md | 2 ++ 2 files changed, 39 insertions(+) create mode 100644 docs/src/fokker_planck_notes.md diff --git a/docs/src/fokker_planck_notes.md b/docs/src/fokker_planck_notes.md new file mode 100644 index 000000000..2be75d955 --- /dev/null +++ b/docs/src/fokker_planck_notes.md @@ -0,0 +1,37 @@ +Fokker Planck collision operator +=============================================== + +We implement the nonlinear Fokker-Planck collision operator for self collisions +using the weak-form finite-element method. This is documented in the +[ExCALIBUR/NEPTUNE report 2070839-TN-07](https://excalibur-neptune.github.io/Documents/TN-07_AHigherOrderFiniteElementImplementationFullFlandauFokkerPlanckCollisionOperatorC.html). +A publication based on this report is in progress. Full online documentation will follow. + +Input parameters +=============================================== + +A series of 0D2V Fokker-Planck input files can be found in + + examples/fokker-planck/ + +and examples of 1D2V pre-sheath simulations with the Fokker-Planck collision operator +can be found in + + examples/fokker-planck-1D2V/ + +noting that the timestepping or resolution parameters may require modification to find +a converged simulation. + +The basic input namelist is structured as follows +``` +[fokker_planck_collisions] +use_fokker_planck = true +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +# for frequency_option = "manual" +nuii = 1.0 +frequency_option = "manual" +``` +Set `use_fokker_planck=false` to turn off Fokker-Planck collisions +without commenting out the namelist. +The default option for `frequency_option = "reference_parameters"`, where `nuii` is set +by the reference parameter inputs. Further specialised input parameters can be +seen in the source at `setup_fkpl_collisions_input()` in `moment_kinetics/src/fokker_planck.jl`. \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 39410aaaf..581c879cc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -13,6 +13,8 @@ Pages = ["getting_started.md", "moment_constraints_notes.md", "boundary_conditions_notes.md", "external_sources_notes.md", + "fokker_planck_notes.md", + "geometry.md", "debugging-hints.md", "developing.md", "manual_setup.md",