Propagate multipole option so it can be used in time evolving code.

moment_kinetics/src/fokker_planck.jl

@@ -94,6 +94,7 @@ function setup_fkpl_collisions_input(toml_input::Dict)
        frequency_option = "reference_parameters",
        self_collisions = true,
        use_conserving_corrections = true,
+       multipole_boundary_data = false,
        slowing_down_test = false,
        sd_density = 1.0,
        sd_temp = 0.01,
@@ -337,13 +338,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
     # 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)
+            vperp, vpa, vperp_spectral, vpa_spectral, 
+            multipole_boundary_data = multipole_boundary_data)
         # 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
@@ -398,7 +401,8 @@ with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi
                          use_Maxwellian_Rosenbluth_coefficients=false,
                          use_Maxwellian_field_particle_distribution=false,
                          algebraic_solve_for_d2Gdvperp2 = false, calculate_GG=false,
-                         calculate_dGdvperp=false) = begin
+                         calculate_dGdvperp=false,
+                         multipole_boundary_data=false) = 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))
@@ -448,7 +452,8 @@ with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi
              d2Gdvpa2,dGdvperp,d2Gdvperpdvpa,d2Gdvperp2,ffsp_in,
              vpa,vperp,vpa_spectral,vperp_spectral,fkpl_arrays,
              algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2,
-             calculate_GG=calculate_GG,calculate_dGdvperp=calculate_dGdvperp)
+             calculate_GG=calculate_GG,calculate_dGdvperp=calculate_dGdvperp,
+             multipole_boundary_data=multipole_boundary_data)
     end
     # assemble the RHS of the collision operator matrix eq
     if use_Maxwellian_field_particle_distribution

moment_kinetics/src/fokker_planck_calculus.jl

@@ -2914,7 +2914,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=false)
+             multipole_boundary_data=false)
 
     # extract the necessary precalculated and buffer arrays from fokkerplanck_arrays
     MM2D_sparse = fkpl_arrays.MM2D_sparse
@@ -2943,7 +2943,7 @@ function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvpe
     rhsvpavperp_copy3 = fkpl_arrays.rhsvpavperp_copy3
 
     # calculate the boundary data
-    if multipole
+    if multipole_boundary_data
         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

@@ -489,6 +489,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
     # 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,9 @@ 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
-        fp_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral; precompute_weights=true)
+        precompute_weights = true && !(collisions.fkpl.multipole_boundary_data)
+        fp_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral;
+                      precompute_weights=precompute_weights)
     else
         fp_arrays = nothing
     end

test_scripts/2D_FEM_assembly_test.jl

@@ -126,9 +126,9 @@ end
         nc_global = vpa.n*vperp.n
         begin_serial_region()
         start_init_time = now()
-
+        precompute_weights = true && !(use_multipole)
         fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral; 
-                           precompute_weights=true, test_dense_matrix_construction=test_dense_construction)
+                           precompute_weights=precompute_weights, test_dense_matrix_construction=test_dense_construction)
         KKpar2D_with_BC_terms_sparse = fkpl_arrays.KKpar2D_with_BC_terms_sparse
         KKperp2D_with_BC_terms_sparse = fkpl_arrays.KKperp2D_with_BC_terms_sparse
         lu_obj_MM = fkpl_arrays.lu_obj_MM
@@ -264,7 +264,8 @@ end
                                              use_Maxwellian_Rosenbluth_coefficients=use_Maxwellian_Rosenbluth_coefficients,
                                              use_Maxwellian_field_particle_distribution=use_Maxwellian_field_particle_distribution,
                                              algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2,
-                                             calculate_GG = false, calculate_dGdvperp=false)
+                                             calculate_GG = false, calculate_dGdvperp=false,
+                                             multipole_boundary_data=use_multipole)
         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)
@@ -275,7 +276,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=use_multipole)
+             algebraic_solve_for_d2Gdvperp2=false,calculate_GG=true,calculate_dGdvperp=true,multipole_boundary_data=use_multipole)
         # extract C[Fs,Fs'] result
         # and Rosenbluth potentials for testing
         begin_s_r_z_anyv_region()