Merge pull request #289 from mabarnes/improve-fkpl-docs

Improvement of online docs for Fokker-Planck options

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`.

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",

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. 
 
@@ -67,7 +67,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
@@ -78,10 +78,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::AbstractDict)
     reference_params = setup_reference_parameters(toml_input)
@@ -153,7 +153,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.
@@ -239,7 +239,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,
@@ -306,7 +307,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,
@@ -504,11 +505,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,
@@ -591,11 +591,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
@@ -616,11 +623,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)
@@ -641,6 +654,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
@@ -688,6 +713,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
@@ -735,7 +769,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
@@ -785,7 +819,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