diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 9a19affb9..a31d3f78e 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/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 diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index b76d553fa..79f6e2ede 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/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 diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index f427c3a83..8d421fa48 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/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 diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index b9d85420d..31a4777a4 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/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 diff --git a/test_scripts/2D_FEM_assembly_test.jl b/test_scripts/2D_FEM_assembly_test.jl index 53c7f9f8e..b6f8408fb 100644 --- a/test_scripts/2D_FEM_assembly_test.jl +++ b/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()