Add variant of multipole_expansion option delta_f_multipole, to use t…

…he exact results for the Rosenbluth potentials for the Maxwellian piece of F, and multipole for the rest. Tests extended to this option. Rename boundary_data -> boundary_data_option to make internal code easier to follow.

moment_kinetics/src/fokker_planck.jl

@@ -95,7 +95,7 @@ function setup_fkpl_collisions_input(toml_input::Dict)
        frequency_option = "reference_parameters",
        self_collisions = true,
        use_conserving_corrections = true,
-       boundary_data = direct_integration,
+       boundary_data_option = direct_integration,
        slowing_down_test = false,
        sd_density = 1.0,
        sd_temp = 0.01,
@@ -339,7 +339,7 @@ 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
-    boundary_data_option = collisions.fkpl.boundary_data
+    boundary_data_option = collisions.fkpl.boundary_data_option
     # N.B. parallelisation using special 'anyv' region
     begin_s_r_z_anyv_region()
     @loop_s_r_z is ir iz begin

moment_kinetics/src/fokker_planck_calculus.jl

@@ -18,6 +18,7 @@ export assemble_explicit_collision_operator_rhs_parallel_analytical_inputs!
 export YY_collision_operator_arrays, calculate_YY_arrays
 export calculate_rosenbluth_potential_boundary_data!
 export calculate_rosenbluth_potential_boundary_data_multipole!
+export calculate_rosenbluth_potential_boundary_data_delta_f_multipole!
 export elliptic_solve!, algebraic_solve!
 export fokkerplanck_arrays_direct_integration_struct
 export fokkerplanck_weakform_arrays_struct
@@ -42,7 +43,10 @@ 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 ..velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure
+using ..input_structs: direct_integration, multipole_expansion, delta_f_multipole
+using ..fokker_planck_test: F_Maxwellian, G_Maxwellian, H_Maxwellian, dHdvpa_Maxwellian, dHdvperp_Maxwellian
+using ..fokker_planck_test: d2Gdvpa2_Maxwellian, d2Gdvperp2_Maxwellian, d2Gdvperpdvpa_Maxwellian, dGdvperp_Maxwellian
 using moment_kinetics.gauss_legendre: get_QQ_local!
 using Dates
 using SpecialFunctions: ellipk, ellipe
@@ -1779,6 +1783,88 @@ function calculate_rosenbluth_potential_boundary_data_multipole!(rpbd::rosenblut
     return nothing
 end
 
+"""
+Function to use the multipole expansion of the Rosenbluth potentials to calculate and
+assign boundary data to an instance of `rosenbluth_potential_boundary_data`, in place,
+without allocation. Use the exact results for the part of F that can be described with 
+a Maxwellian, and the multipole expansion for the remainder.
+"""
+function calculate_rosenbluth_potential_boundary_data_delta_f_multipole!(rpbd::rosenbluth_potential_boundary_data,
+    pdf,dummy_vpavperp,vpa,vperp,vpa_spectral,vperp_spectral;
+    calculate_GG=false,calculate_dGdvperp=false)
+
+    mass = 1.0
+    # first, compute the moments and delta f
+    begin_anyv_region()
+    @anyv_serial_region begin
+      dens =  get_density(pdf, vpa, vperp)
+      upar = get_upar(pdf, vpa, vperp, dens)
+      ppar = get_ppar(pdf, vpa, vperp, upar)
+      pperp = get_pperp(pdf, vpa, vperp)
+      pressure = get_pressure(ppar,pperp)
+      vth = sqrt(2.0*pressure/(dens*mass))
+      @loop_vperp_vpa ivperp ivpa begin
+          dummy_vpavperp[ivpa,ivperp] = pdf[ivpa,ivperp] - F_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) 
+      end
+    end
+    # ensure data is synchronized
+    _anyv_subblock_synchronize()
+    # now pass the delta f to the multipole function
+    calculate_rosenbluth_potential_boundary_data_multipole!(rpbd,dummy_vpavperp,
+      vpa,vperp,vpa_spectral,vperp_spectral,
+      calculate_GG=calculate_GG,calculate_dGdvperp=calculate_dGdvperp)
+    # now add on the contributions from the Maxwellian
+    nvpa = vpa.n
+    nvperp = vperp.n
+    begin_anyv_vperp_region()
+    @loop_vperp ivperp begin
+                rpbd.H_data.lower_boundary_vpa[ivperp] += H_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                rpbd.H_data.upper_boundary_vpa[ivperp] += H_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)
+                rpbd.dHdvpa_data.lower_boundary_vpa[ivperp] += dHdvpa_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                rpbd.dHdvpa_data.upper_boundary_vpa[ivperp] += dHdvpa_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)
+                rpbd.dHdvperp_data.lower_boundary_vpa[ivperp] += dHdvperp_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                rpbd.dHdvperp_data.upper_boundary_vpa[ivperp] += dHdvperp_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)
+                rpbd.d2Gdvpa2_data.lower_boundary_vpa[ivperp] += d2Gdvpa2_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                rpbd.d2Gdvpa2_data.upper_boundary_vpa[ivperp] += d2Gdvpa2_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)                
+                rpbd.d2Gdvperpdvpa_data.lower_boundary_vpa[ivperp] += d2Gdvperpdvpa_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                rpbd.d2Gdvperpdvpa_data.upper_boundary_vpa[ivperp] += d2Gdvperpdvpa_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)                
+                rpbd.d2Gdvperp2_data.lower_boundary_vpa[ivperp] += d2Gdvperp2_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                rpbd.d2Gdvperp2_data.upper_boundary_vpa[ivperp] += d2Gdvperp2_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)
+    end
+    begin_anyv_vpa_region()
+    @loop_vpa ivpa begin
+                rpbd.H_data.upper_boundary_vperp[ivpa] += H_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+                rpbd.dHdvpa_data.upper_boundary_vperp[ivpa] += dHdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+                rpbd.dHdvperp_data.upper_boundary_vperp[ivpa] += dHdvperp_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+                rpbd.d2Gdvpa2_data.upper_boundary_vperp[ivpa] += d2Gdvpa2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+                rpbd.d2Gdvperpdvpa_data.upper_boundary_vperp[ivpa] += d2Gdvperpdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+                rpbd.d2Gdvperp2_data.upper_boundary_vperp[ivpa] += d2Gdvperp2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+    end
+    if calculate_GG
+       begin_anyv_vperp_region()
+       @loop_vperp ivperp begin
+                   rpbd.G_data.lower_boundary_vpa[ivperp] += G_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                   rpbd.G_data.upper_boundary_vpa[ivperp] += G_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)
+       end
+       begin_anyv_vpa_region()
+       @loop_vpa ivpa begin
+                   rpbd.G_data.upper_boundary_vperp[ivpa] += G_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+       end
+    end
+    if calculate_dGdvperp
+       begin_anyv_vperp_region()
+       @loop_vperp ivperp begin
+                   rpbd.dGdvperp_data.lower_boundary_vpa[ivperp] += dGdvperp_Maxwellian(dens,upar,vth,vpa,vperp,1,ivperp)
+                   rpbd.dGdvperp_data.upper_boundary_vpa[ivperp] += dGdvperp_Maxwellian(dens,upar,vth,vpa,vperp,nvpa,ivperp)
+       end
+       begin_anyv_vpa_region()
+       @loop_vpa ivpa begin
+                   rpbd.dGdvperp_data.upper_boundary_vperp[ivpa] += dGdvperp_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,nvperp)
+       end
+    end
+    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
@@ -2834,9 +2920,17 @@ function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvpe
     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
+    elseif boundary_data_option == delta_f_multipole # use a variant of the multipole method
+        calculate_rosenbluth_potential_boundary_data_delta_f_multipole!(rpbd,ffsp_in,S_dummy,vpa,vperp,vpa_spectral,vperp_spectral,
+          calculate_GG=calculate_GG,calculate_dGdvperp=(calculate_dGdvperp||algebraic_solve_for_d2Gdvperp2))
+    elseif boundary_data_option == direct_integration  # use direct integration on the boundary
         calculate_rosenbluth_potential_boundary_data!(rpbd,bwgt,ffsp_in,vpa,vperp,vpa_spectral,vperp_spectral,
          calculate_GG=calculate_GG,calculate_dGdvperp=(calculate_dGdvperp||algebraic_solve_for_d2Gdvperp2))
+    else
+        error("No valid boundary_data_option specified. \n 
+              Pick  boundary_data_option='$multipole_expansion' \n 
+              or  boundary_data_option='$delta_f_multipole' \n 
+              or boundary_data_option='$direct_integration'")
     end
     # carry out the elliptic solves required
     begin_anyv_vperp_vpa_region()

moment_kinetics/src/input_structs.jl

@@ -482,10 +482,12 @@ end
 @enum boundary_data_type begin
     direct_integration
     multipole_expansion
+    delta_f_multipole
 end
 export boundary_data_type
 export direct_integration
 export multipole_expansion
+export delta_f_multipole
 
 Base.@kwdef struct fkpl_collisions_input
     # option to check if fokker planck frequency should be > 0
@@ -500,7 +502,7 @@ Base.@kwdef struct fkpl_collisions_input
     # option to determine if ad-hoc moment_kinetics-style conserving corrections are used
     use_conserving_corrections::Bool
     # enum option to determine which method is used to provide boundary data for Rosenbluth potential calculations.
-    boundary_data::boundary_data_type
+    boundary_data_option::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,7 @@ 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
-        if collisions.fkpl.boundary_data == direct_integration
+        if collisions.fkpl.boundary_data_option == direct_integration
             precompute_weights = true
         else
             precompute_weights = false

moment_kinetics/test/fokker_planck_tests.jl

@@ -11,7 +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.input_structs: direct_integration, multipole_expansion, delta_f_multipole
 
 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
@@ -208,7 +208,7 @@ function runtests()
 
         @testset "weak-form Rosenbluth potential calculation: elliptic solve" begin
             println("    - test weak-form Rosenbluth potential calculation: elliptic solve")
-            @testset "$boundary_data_option" for boundary_data_option in (direct_integration,multipole_expansion)
+            @testset "$boundary_data_option" for boundary_data_option in (direct_integration,multipole_expansion,delta_f_multipole)
                 println("        -  boundary_data_option=$boundary_data_option")
                 ngrid = 9
                 nelement_vpa = 8