Merge pull request #291 from mabarnes/multipole-boundary-data

Multipole boundary data

moment_kinetics/src/external_sources.jl

@@ -814,7 +814,7 @@ Add external source term to the ion kinetic equation.
         end
     elseif source_type == "alphas" || source_type == "alphas-with-losses"
         begin_s_r_z_region()
-        source_v0 = ion_source_settings.source_v0
+        source_v0 = ion_source.source_v0
         if !(source_v0 > 1.0e-8)
             error("source_v0=$source_v0 < 1.0e-8")
         end
@@ -840,8 +840,8 @@ Add external source term to the ion kinetic equation.
             end
 
             if source_type == "alphas-with-losses"
-                sink_vth = ion_source_settings.sink_vth
-                sink_strength = ion_source_settings.sink_strength
+                sink_vth = ion_source.sink_vth
+                sink_strength = ion_source.sink_strength
                 if !(sink_vth > 1.0e-8)
                    error("sink_vth=$sink_vth < 1.0e-8")
                 end
@@ -871,8 +871,8 @@ Add external source term to the ion kinetic equation.
         end
     elseif source_type == "beam" || source_type == "beam-with-losses"
         begin_s_r_z_region()
-        source_vpa0 = ion_source_settings.source_vpa0
-        source_vperp0 = ion_source_settings.source_vperp0
+        source_vpa0 = ion_source.source_vpa0
+        source_vperp0 = ion_source.source_vperp0
         if !(source_vpa0 > 1.0e-8)
             error("source_vpa0=$source_vpa0 < 1.0e-8")
         end
@@ -901,8 +901,8 @@ Add external source term to the ion kinetic equation.
             end
 
             if source_type == "beam-with-losses"
-                sink_vth = ion_source_settings.sink_vth
-                sink_strength = ion_source_settings.sink_strength
+                sink_vth = ion_source.sink_vth
+                sink_strength = ion_source.sink_strength
                 if !(sink_vth > 1.0e-8)
                    error("sink_vth=$sink_vth < 1.0e-8")
                 end

moment_kinetics/src/fokker_planck.jl

@@ -54,6 +54,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!
@@ -95,6 +96,7 @@ function setup_fkpl_collisions_input(toml_input::AbstractDict)
        frequency_option = "reference_parameters",
        self_collisions = true,
        use_conserving_corrections = true,
+       boundary_data_option = direct_integration,
        slowing_down_test = false,
        sd_density = 1.0,
        sd_temp = 0.01,
@@ -338,13 +340,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
+    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
         # 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, 
+            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
@@ -399,7 +403,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,
+                         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))
@@ -449,7 +454,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,
+             boundary_data_option=boundary_data_option)
     end
     # assemble the RHS of the collision operator matrix eq
     if use_Maxwellian_field_particle_distribution