Add a more descriptive enum option to switch between boundary data op…

…tions for the calculation of the Rosenbluth potentials.

moment_kinetics/src/fokker_planck.jl

@@ -53,6 +53,7 @@ using ..velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_qpar,
 using ..looping
 using ..timer_utils
 using ..input_structs: fkpl_collisions_input, set_defaults_and_check_section!
+using ..input_structs: multipole_expansion, direct_integration
 using ..reference_parameters: get_reference_collision_frequency_ii
 using ..fokker_planck_calculus: init_Rosenbluth_potential_integration_weights!
 using ..fokker_planck_calculus: init_Rosenbluth_potential_boundary_integration_weights!
@@ -94,7 +95,7 @@ function setup_fkpl_collisions_input(toml_input::Dict)
        frequency_option = "reference_parameters",
        self_collisions = true,
        use_conserving_corrections = true,
-       multipole_boundary_data = false,
+       boundary_data = direct_integration,
        slowing_down_test = false,
        sd_density = 1.0,
        sd_temp = 0.01,
@@ -338,15 +339,15 @@ Function for advancing with the explicit, weak-form, self-collision operator.
     Zi = collisions.fkpl.Zi # generalise!
     nussp = nuref*(Zi^4) # include charge number factor for self collisions
     use_conserving_corrections = collisions.fkpl.use_conserving_corrections
-    multipole_boundary_data = collisions.fkpl.multipole_boundary_data
+    boundary_data_option = collisions.fkpl.boundary_data
     # N.B. parallelisation using special 'anyv' region
     begin_s_r_z_anyv_region()
     @loop_s_r_z is ir iz begin
         # first argument is Fs, and second argument is Fs' in C[Fs,Fs'] 
         @views fokker_planck_collision_operator_weak_form!(
             pdf_in[:,:,iz,ir,is], pdf_in[:,:,iz,ir,is], ms, msp, nussp, fkpl_arrays,
             vperp, vpa, vperp_spectral, vpa_spectral, 
-            multipole_boundary_data = multipole_boundary_data)
+            boundary_data_option = boundary_data_option)
         # enforce the boundary conditions on CC before it is used for timestepping
         enforce_vpavperp_BCs!(fkpl_arrays.CC,vpa,vperp,vpa_spectral,vperp_spectral)
         # make ad-hoc conserving corrections
@@ -402,7 +403,7 @@ with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi
                          use_Maxwellian_field_particle_distribution=false,
                          algebraic_solve_for_d2Gdvperp2 = false, calculate_GG=false,
                          calculate_dGdvperp=false,
-                         multipole_boundary_data=false) = begin
+                         boundary_data_option=direct_integration) = begin
     @boundscheck vpa.n == size(ffsp_in,1) || throw(BoundsError(ffsp_in))
     @boundscheck vperp.n == size(ffsp_in,2) || throw(BoundsError(ffsp_in))
     @boundscheck vpa.n == size(ffs_in,1) || throw(BoundsError(ffs_in))
@@ -453,7 +454,7 @@ with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi
              vpa,vperp,vpa_spectral,vperp_spectral,fkpl_arrays,
              algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2,
              calculate_GG=calculate_GG,calculate_dGdvperp=calculate_dGdvperp,
-             multipole_boundary_data=multipole_boundary_data)
+             boundary_data_option=boundary_data_option)
     end
     # assemble the RHS of the collision operator matrix eq
     if use_Maxwellian_field_particle_distribution

moment_kinetics/src/fokker_planck_calculus.jl

@@ -42,6 +42,7 @@ using ..communication: MPISharedArray, global_rank
 using ..lagrange_polynomials: lagrange_poly, lagrange_poly_optimised
 using ..looping
 using ..velocity_moments: integrate_over_vspace
+using ..input_structs: direct_integration, multipole_expansion
 using moment_kinetics.gauss_legendre: get_QQ_local!
 using Dates
 using SpecialFunctions: ellipk, ellipe
@@ -2914,7 +2915,7 @@ function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvpe
              d2Gdvpa2,dGdvperp,d2Gdvperpdvpa,d2Gdvperp2,ffsp_in,
              vpa,vperp,vpa_spectral,vperp_spectral,fkpl_arrays::fokkerplanck_weakform_arrays_struct;
              algebraic_solve_for_d2Gdvperp2=false,calculate_GG=false,calculate_dGdvperp=false,
-             multipole_boundary_data=false)
+             boundary_data_option=direct_integration)
 
     # extract the necessary precalculated and buffer arrays from fokkerplanck_arrays
     MM2D_sparse = fkpl_arrays.MM2D_sparse
@@ -2943,7 +2944,7 @@ function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvpe
     rhsvpavperp_copy3 = fkpl_arrays.rhsvpavperp_copy3
 
     # calculate the boundary data
-    if multipole_boundary_data
+    if boundary_data_option == multipole_expansion
         calculate_rosenbluth_potential_boundary_data_multipole!(rpbd,ffsp_in,vpa,vperp,vpa_spectral,vperp_spectral,
           calculate_GG=calculate_GG,calculate_dGdvperp=(calculate_dGdvperp||algebraic_solve_for_d2Gdvperp2))
     else # use direct integration on the boundary

moment_kinetics/src/input_structs.jl

@@ -477,6 +477,16 @@ Base.@kwdef struct krook_collisions_input
     frequency_option::String # "reference_parameters" # "manual", 
 end
 
+"""
+"""
+@enum boundary_data_type begin
+    direct_integration
+    multipole_expansion
+end
+export boundary_data_type
+export direct_integration
+export multipole_expansion
+
 Base.@kwdef struct fkpl_collisions_input
     # option to check if fokker planck frequency should be > 0
     use_fokker_planck::Bool
@@ -489,8 +499,8 @@ Base.@kwdef struct fkpl_collisions_input
     self_collisions::Bool
     # option to determine if ad-hoc moment_kinetics-style conserving corrections are used
     use_conserving_corrections::Bool
-    # option to determine if multipole expansion is used to provide boundary data for Rosenbluth potential calculations.
-    multipole_boundary_data::Bool
+    # enum option to determine which method is used to provide boundary data for Rosenbluth potential calculations.
+    boundary_data::boundary_data_type
     # option to determine if cross-collisions against fixed Maxwellians are used
     slowing_down_test::Bool
     # Setting to switch between different options for Fokker-Planck collision frequency input

moment_kinetics/src/time_advance.jl

@@ -734,7 +734,11 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop
                                                   n_neutral_species_alloc, t_params)
     # create arrays for Fokker-Planck collisions 
     if advance.explicit_weakform_fp_collisions
-        precompute_weights = true && !(collisions.fkpl.multipole_boundary_data)
+        if collisions.fkpl.boundary_data == direct_integration
+            precompute_weights = true
+        else
+            precompute_weights = false
+        end
         fp_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral;
                       precompute_weights=precompute_weights)
     else

moment_kinetics/test/fokker_planck_tests.jl

@@ -11,6 +11,7 @@ using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
 using moment_kinetics.coordinates: define_coordinate
 using moment_kinetics.type_definitions: mk_float, mk_int
 using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure
+using moment_kinetics.input_structs: direct_integration, multipole_expansion
 
 using moment_kinetics.fokker_planck: init_fokker_planck_collisions_weak_form, fokker_planck_collision_operator_weak_form!
 using moment_kinetics.fokker_planck: conserving_corrections!, init_fokker_planck_collisions_direct_integration
@@ -207,16 +208,21 @@ function runtests()
 
         @testset "weak-form Rosenbluth potential calculation: elliptic solve" begin
             println("    - test weak-form Rosenbluth potential calculation: elliptic solve")
-            @testset "$multipole_boundary_data" for multipole_boundary_data in (true,false)
-                println("        -  multipole_boundary_data=$multipole_boundary_data")
+            @testset "$boundary_data_option" for boundary_data_option in (direct_integration,multipole_expansion)
+                println("        -  boundary_data_option=$boundary_data_option")
                 ngrid = 9
                 nelement_vpa = 8
                 nelement_vperp = 4
                 vpa, vpa_spectral, vperp, vperp_spectral = create_grids(ngrid,nelement_vpa,nelement_vperp,
                                                                             Lvpa=12.0,Lvperp=6.0)
                 begin_serial_region()
+                if boundary_data_option == direct_integration
+                    precompute_weights = true
+                else
+                    precompute_weights = false
+                end
                 fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral,
-                                                                      precompute_weights=(true &&!(multipole_boundary_data)),
+                                                                      precompute_weights=precompute_weights,
                                                                       print_to_screen=print_to_screen)
                 dummy_array = allocate_float(vpa.n,vperp.n)
                 F_M = allocate_float(vpa.n,vperp.n)
@@ -274,7 +280,7 @@ function runtests()
                      fkpl_arrays.d2Gdvperp2, F_M, vpa, vperp, vpa_spectral, vperp_spectral,
                      fkpl_arrays; algebraic_solve_for_d2Gdvperp2=false,
                      calculate_GG=true, calculate_dGdvperp=true,
-                     multipole_boundary_data=multipole_boundary_data)
+                     boundary_data_option=boundary_data_option)
                 # extract C[Fs,Fs'] result
                 # and Rosenbluth potentials for testing
                 begin_s_r_z_anyv_region()
@@ -296,7 +302,7 @@ function runtests()
                     max_dHdvperp_boundary_data_err, max_G_boundary_data_err,
                     max_dGdvperp_boundary_data_err, max_d2Gdvperp2_boundary_data_err,
                     max_d2Gdvperpdvpa_boundary_data_err, max_d2Gdvpa2_boundary_data_err = test_rosenbluth_potential_boundary_data(fkpl_arrays.rpbd,rpbd_exact,vpa,vperp,print_to_screen=print_to_screen)
-                    if multipole_boundary_data
+                    if boundary_data_option==multipole_expansion
                         atol_max_H = 5.0e-8
                         atol_max_dHdvpa = 5.0e-8
                         atol_max_dHdvperp = 5.0e-8
@@ -332,7 +338,7 @@ function runtests()
                     dGdvperp_M_max, dGdvperp_M_L2 = print_test_data(dGdvperp_M_exact,dGdvperp_M_num,dGdvperp_M_err,"dGdvperp_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen)
                     d2Gdvperpdvpa_M_max, d2Gdvperpdvpa_M_L2 = print_test_data(d2Gdvperpdvpa_M_exact,d2Gdvperpdvpa_M_num,d2Gdvperpdvpa_M_err,"d2Gdvperpdvpa_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen)
                     d2Gdvperp2_M_max, d2Gdvperp2_M_L2 = print_test_data(d2Gdvperp2_M_exact,d2Gdvperp2_M_num,d2Gdvperp2_M_err,"d2Gdvperp2_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen)
-                    if multipole_boundary_data
+                    if boundary_data_option==multipole_expansion
                         atol_max_H = 2.0e-7
                         atol_L2_H = 5.0e-9
                         atol_max_dHdvpa = 2.0e-6

test_scripts/2D_FEM_assembly_test.jl

@@ -19,6 +19,7 @@ using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_ppe
 using moment_kinetics.communication
 using moment_kinetics.communication: MPISharedArray
 using moment_kinetics.looping
+using moment_kinetics.input_structs: direct_integration, multipole_expansion
 using SparseArrays: sparse
 using LinearAlgebra: mul!, lu, cholesky
 
@@ -79,7 +80,7 @@ end
         use_Maxwellian_field_particle_distribution=false,
         test_numerical_conserving_terms=false,
         algebraic_solve_for_d2Gdvperp2=false,
-        use_multipole=false)
+        boundary_data_option=direct_integration)
         # define inputs needed for the test
         #plot_test_output = false#true
         #test_parallelism = false#true
@@ -126,7 +127,11 @@ end
         nc_global = vpa.n*vperp.n
         begin_serial_region()
         start_init_time = now()
-        precompute_weights = true && !(use_multipole)
+        if boundary_data_option == direct_integration
+            precompute_weights = true
+        else
+            precompute_weights = false
+        end
         fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral; 
                            precompute_weights=precompute_weights, test_dense_matrix_construction=test_dense_construction)
         KKpar2D_with_BC_terms_sparse = fkpl_arrays.KKpar2D_with_BC_terms_sparse
@@ -265,7 +270,7 @@ end
                                              use_Maxwellian_field_particle_distribution=use_Maxwellian_field_particle_distribution,
                                              algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2,
                                              calculate_GG = false, calculate_dGdvperp=false,
-                                             multipole_boundary_data=use_multipole)
+                                             boundary_data_option=boundary_data_option)
         if test_numerical_conserving_terms && test_self_operator
             # enforce the boundary conditions on CC before it is used for timestepping
             enforce_vpavperp_BCs!(fkpl_arrays.CC,vpa,vperp,vpa_spectral,vperp_spectral)
@@ -276,7 +281,7 @@ end
         calculate_rosenbluth_potentials_via_elliptic_solve!(fkpl_arrays.GG,fkpl_arrays.HH,fkpl_arrays.dHdvpa,fkpl_arrays.dHdvperp,
              fkpl_arrays.d2Gdvpa2,fkpl_arrays.dGdvperp,fkpl_arrays.d2Gdvperpdvpa,fkpl_arrays.d2Gdvperp2,F_M,
              vpa,vperp,vpa_spectral,vperp_spectral,fkpl_arrays;
-             algebraic_solve_for_d2Gdvperp2=false,calculate_GG=true,calculate_dGdvperp=true,multipole_boundary_data=use_multipole)
+             algebraic_solve_for_d2Gdvperp2=false,calculate_GG=true,calculate_dGdvperp=true,boundary_data_option=boundary_data_option)
         # extract C[Fs,Fs'] result
         # and Rosenbluth potentials for testing
         begin_s_r_z_anyv_region()
@@ -392,7 +397,7 @@ end
         algebraic_solve_for_d2Gdvperp2=false,
         test_self_operator = true,
         Lvpa = 12.0, Lvperp = 6.0,
-        use_multipole = false)
+        boundary_data_option = direct_integration)
         initialize_comms!()
         #ngrid = 5
         #plot_scan = true
@@ -463,7 +468,7 @@ end
             use_Maxwellian_field_particle_distribution=use_Maxwellian_field_particle_distribution,
             test_numerical_conserving_terms=test_numerical_conserving_terms,
             algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2,
-            standalone=false, Lvpa=Lvpa, Lvperp=Lvperp, use_multipole=use_multipole)
+            standalone=false, Lvpa=Lvpa, Lvperp=Lvperp, boundary_data_option=boundary_data_option)
             max_C_err[iscan], L2_C_err[iscan] = fkerr.C_M.max ,fkerr.C_M.L2
             max_H_err[iscan], L2_H_err[iscan] = fkerr.H_M.max ,fkerr.H_M.L2
             max_dHdvpa_err[iscan], L2_dHdvpa_err[iscan] = fkerr.dHdvpa_M.max ,fkerr.dHdvpa_M.L2