From 59fd6949279388e3f34f48097d00ac5eec0e1e10 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Sat, 9 Nov 2024 20:40:34 +0000 Subject: [PATCH] Add a more descriptive enum option to switch between boundary data options for the calculation of the Rosenbluth potentials. --- moment_kinetics/src/fokker_planck.jl | 11 ++++++----- moment_kinetics/src/fokker_planck_calculus.jl | 5 +++-- moment_kinetics/src/input_structs.jl | 14 ++++++++++++-- moment_kinetics/src/time_advance.jl | 6 +++++- moment_kinetics/test/fokker_planck_tests.jl | 18 ++++++++++++------ test_scripts/2D_FEM_assembly_test.jl | 17 +++++++++++------ 6 files changed, 49 insertions(+), 22 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index a31d3f78e..168a222b3 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -53,6 +53,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! @@ -94,7 +95,7 @@ function setup_fkpl_collisions_input(toml_input::Dict) frequency_option = "reference_parameters", self_collisions = true, use_conserving_corrections = true, - multipole_boundary_data = false, + boundary_data = direct_integration, slowing_down_test = false, sd_density = 1.0, sd_temp = 0.01, @@ -338,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 - multipole_boundary_data = collisions.fkpl.multipole_boundary_data + boundary_data_option = collisions.fkpl.boundary_data # N.B. parallelisation using special 'anyv' region begin_s_r_z_anyv_region() @loop_s_r_z is ir iz begin @@ -346,7 +347,7 @@ Function for advancing with the explicit, weak-form, self-collision operator. @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, - multipole_boundary_data = multipole_boundary_data) + 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 @@ -402,7 +403,7 @@ with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi use_Maxwellian_field_particle_distribution=false, algebraic_solve_for_d2Gdvperp2 = false, calculate_GG=false, calculate_dGdvperp=false, - multipole_boundary_data=false) = begin + 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)) @@ -453,7 +454,7 @@ with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi vpa,vperp,vpa_spectral,vperp_spectral,fkpl_arrays, algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2, calculate_GG=calculate_GG,calculate_dGdvperp=calculate_dGdvperp, - multipole_boundary_data=multipole_boundary_data) + boundary_data_option=boundary_data_option) 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 79f6e2ede..9831927f4 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -42,6 +42,7 @@ 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 moment_kinetics.gauss_legendre: get_QQ_local! using Dates using SpecialFunctions: ellipk, ellipe @@ -2914,7 +2915,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_boundary_data=false) + boundary_data_option=direct_integration) # extract the necessary precalculated and buffer arrays from fokkerplanck_arrays MM2D_sparse = fkpl_arrays.MM2D_sparse @@ -2943,7 +2944,7 @@ function calculate_rosenbluth_potentials_via_elliptic_solve!(GG,HH,dHdvpa,dHdvpe rhsvpavperp_copy3 = fkpl_arrays.rhsvpavperp_copy3 # calculate the boundary data - if multipole_boundary_data + 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 diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 8d421fa48..8a0fe2388 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -477,6 +477,16 @@ Base.@kwdef struct krook_collisions_input frequency_option::String # "reference_parameters" # "manual", end +""" +""" +@enum boundary_data_type begin + direct_integration + multipole_expansion +end +export boundary_data_type +export direct_integration +export multipole_expansion + Base.@kwdef struct fkpl_collisions_input # option to check if fokker planck frequency should be > 0 use_fokker_planck::Bool @@ -489,8 +499,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 + # enum option to determine which method is used to provide boundary data for Rosenbluth potential calculations. + boundary_data::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 31a4777a4..3d4ec26e8 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -734,7 +734,11 @@ 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 - precompute_weights = true && !(collisions.fkpl.multipole_boundary_data) + if collisions.fkpl.boundary_data == direct_integration + precompute_weights = true + else + precompute_weights = false + end fp_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral; precompute_weights=precompute_weights) else diff --git a/moment_kinetics/test/fokker_planck_tests.jl b/moment_kinetics/test/fokker_planck_tests.jl index 16fad7bb0..93b6d687e 100644 --- a/moment_kinetics/test/fokker_planck_tests.jl +++ b/moment_kinetics/test/fokker_planck_tests.jl @@ -11,6 +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.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 @@ -207,16 +208,21 @@ function runtests() @testset "weak-form Rosenbluth potential calculation: elliptic solve" begin println(" - test weak-form Rosenbluth potential calculation: elliptic solve") - @testset "$multipole_boundary_data" for multipole_boundary_data in (true,false) - println(" - multipole_boundary_data=$multipole_boundary_data") + @testset "$boundary_data_option" for boundary_data_option in (direct_integration,multipole_expansion) + println(" - boundary_data_option=$boundary_data_option") ngrid = 9 nelement_vpa = 8 nelement_vperp = 4 vpa, vpa_spectral, vperp, vperp_spectral = create_grids(ngrid,nelement_vpa,nelement_vperp, Lvpa=12.0,Lvperp=6.0) begin_serial_region() + if boundary_data_option == direct_integration + precompute_weights = true + else + precompute_weights = false + end fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral, - precompute_weights=(true &&!(multipole_boundary_data)), + precompute_weights=precompute_weights, print_to_screen=print_to_screen) dummy_array = allocate_float(vpa.n,vperp.n) F_M = allocate_float(vpa.n,vperp.n) @@ -274,7 +280,7 @@ function runtests() 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_boundary_data=multipole_boundary_data) + boundary_data_option=boundary_data_option) # extract C[Fs,Fs'] result # and Rosenbluth potentials for testing begin_s_r_z_anyv_region() @@ -296,7 +302,7 @@ function runtests() max_dHdvperp_boundary_data_err, max_G_boundary_data_err, max_dGdvperp_boundary_data_err, max_d2Gdvperp2_boundary_data_err, max_d2Gdvperpdvpa_boundary_data_err, max_d2Gdvpa2_boundary_data_err = test_rosenbluth_potential_boundary_data(fkpl_arrays.rpbd,rpbd_exact,vpa,vperp,print_to_screen=print_to_screen) - if multipole_boundary_data + if boundary_data_option==multipole_expansion atol_max_H = 5.0e-8 atol_max_dHdvpa = 5.0e-8 atol_max_dHdvperp = 5.0e-8 @@ -332,7 +338,7 @@ function runtests() dGdvperp_M_max, dGdvperp_M_L2 = print_test_data(dGdvperp_M_exact,dGdvperp_M_num,dGdvperp_M_err,"dGdvperp_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen) d2Gdvperpdvpa_M_max, d2Gdvperpdvpa_M_L2 = print_test_data(d2Gdvperpdvpa_M_exact,d2Gdvperpdvpa_M_num,d2Gdvperpdvpa_M_err,"d2Gdvperpdvpa_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen) d2Gdvperp2_M_max, d2Gdvperp2_M_L2 = print_test_data(d2Gdvperp2_M_exact,d2Gdvperp2_M_num,d2Gdvperp2_M_err,"d2Gdvperp2_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen) - if multipole_boundary_data + if boundary_data_option==multipole_expansion atol_max_H = 2.0e-7 atol_L2_H = 5.0e-9 atol_max_dHdvpa = 2.0e-6 diff --git a/test_scripts/2D_FEM_assembly_test.jl b/test_scripts/2D_FEM_assembly_test.jl index b6f8408fb..69742fb09 100644 --- a/test_scripts/2D_FEM_assembly_test.jl +++ b/test_scripts/2D_FEM_assembly_test.jl @@ -19,6 +19,7 @@ using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_ppe using moment_kinetics.communication using moment_kinetics.communication: MPISharedArray using moment_kinetics.looping +using moment_kinetics.input_structs: direct_integration, multipole_expansion using SparseArrays: sparse using LinearAlgebra: mul!, lu, cholesky @@ -79,7 +80,7 @@ end use_Maxwellian_field_particle_distribution=false, test_numerical_conserving_terms=false, algebraic_solve_for_d2Gdvperp2=false, - use_multipole=false) + boundary_data_option=direct_integration) # define inputs needed for the test #plot_test_output = false#true #test_parallelism = false#true @@ -126,7 +127,11 @@ end nc_global = vpa.n*vperp.n begin_serial_region() start_init_time = now() - precompute_weights = true && !(use_multipole) + if boundary_data_option == direct_integration + precompute_weights = true + else + precompute_weights = false + end fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral; precompute_weights=precompute_weights, test_dense_matrix_construction=test_dense_construction) KKpar2D_with_BC_terms_sparse = fkpl_arrays.KKpar2D_with_BC_terms_sparse @@ -265,7 +270,7 @@ end use_Maxwellian_field_particle_distribution=use_Maxwellian_field_particle_distribution, algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2, calculate_GG = false, calculate_dGdvperp=false, - multipole_boundary_data=use_multipole) + boundary_data_option=boundary_data_option) 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) @@ -276,7 +281,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_boundary_data=use_multipole) + algebraic_solve_for_d2Gdvperp2=false,calculate_GG=true,calculate_dGdvperp=true,boundary_data_option=boundary_data_option) # extract C[Fs,Fs'] result # and Rosenbluth potentials for testing begin_s_r_z_anyv_region() @@ -392,7 +397,7 @@ end algebraic_solve_for_d2Gdvperp2=false, test_self_operator = true, Lvpa = 12.0, Lvperp = 6.0, - use_multipole = false) + boundary_data_option = direct_integration) initialize_comms!() #ngrid = 5 #plot_scan = true @@ -463,7 +468,7 @@ end use_Maxwellian_field_particle_distribution=use_Maxwellian_field_particle_distribution, test_numerical_conserving_terms=test_numerical_conserving_terms, algebraic_solve_for_d2Gdvperp2=algebraic_solve_for_d2Gdvperp2, - standalone=false, Lvpa=Lvpa, Lvperp=Lvperp, use_multipole=use_multipole) + standalone=false, Lvpa=Lvpa, Lvperp=Lvperp, boundary_data_option=boundary_data_option) max_C_err[iscan], L2_C_err[iscan] = fkerr.C_M.max ,fkerr.C_M.L2 max_H_err[iscan], L2_H_err[iscan] = fkerr.H_M.max ,fkerr.H_M.L2 max_dHdvpa_err[iscan], L2_dHdvpa_err[iscan] = fkerr.dHdvpa_M.max ,fkerr.dHdvpa_M.L2