diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 168a222b3..6ac0e0543 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/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 diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index ef78ef293..5594c8e45 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/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() diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 8a0fe2388..84d78bd2c 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/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 diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 3d4ec26e8..f2df9030a 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/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 diff --git a/moment_kinetics/test/fokker_planck_tests.jl b/moment_kinetics/test/fokker_planck_tests.jl index 93b6d687e..317798a05 100644 --- a/moment_kinetics/test/fokker_planck_tests.jl +++ b/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