diff --git a/examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml b/examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml index 757712013..3d4e041e3 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-report-resolutions.toml @@ -30,13 +30,6 @@ charge_exchange_frequency = 0.0 ionization_frequency = 0.0 constant_ionization_rate = false -nstep = 200000 -dt = 1.0e-3 -nwrite = 50 -nwrite_dfns = 50 -use_semi_lagrange = false -n_rk_stages = 4 -split_operators = false z_ngrid = 1 z_nelement = 1 z_nelement_local = 1 @@ -64,3 +57,9 @@ use_fokker_planck = true # nuii sets the normalised input C[F,F] Fokker-Planck collision frequency nuii = 1.0 frequency_option = "manual" + +[timestepping] +nstep = 200000 +dt = 1.0e-3 +nwrite = 50 +nwrite_dfns = 50 diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml new file mode 100644 index 000000000..098dc787a --- /dev/null +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml @@ -0,0 +1,91 @@ +# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons. +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +T_e = 1.0 +T_wall = 1.0 +initial_density1 = 1.0 +initial_temperature1 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.0 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +vpa_IC_option1 = "isotropic-beam" +vpa_IC_v01 = 1.0 +vpa_IC_vth01 = 0.1 +#vpa_IC_option1 = "directed-beam" +#vpa_IC_vpa01 = -1.5 +#vpa_IC_vperp01 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false + +z_ngrid = 1 +z_nelement = 1 +z_nelement_local = 1 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 3.0 +vpa_bc = "zero" +vpa_discretization = "gausslegendre_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 1.5 +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero" +# Fokker-Planck operator requires the "gausslegendre_pseudospectral +# options for the vpa and vperp grids + +[fokker_planck_collisions] +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic) +Zi = 2.0 +self_collisions = true +slowing_down_test = true +frequency_option = "manual" +use_fokker_planck = true +sd_density = 1.0 +sd_temp = 0.0025 # TD/Ealpha +sd_mi = 0.5 # mD/malpha +sd_me = 0.000013616 # 0.25/1836.0 me/malpha +sd_q = 1.0 + +[ion_source] +active=true +source_strength=1.0 +source_T=0.005 +source_n=1.0 +r_profile="constant" +r_width=1.0 +r_relative_minimum=0.0 +z_profile="gaussian" +z_width=0.1 +z_relative_minimum=0.0 +source_v0 = 1.0 +#source_type="alphas" +source_type="alphas-with-losses" +#source_type="beam-with-losses" +#source_vpa0 = 1.0 +#source_vperp0 = 1.0 +sink_strength = 1.0 +sink_vth=0.07071067811865475 + +[timestepping] +nstep = 50000 +dt = 5.0e-4 +nwrite = 100 +nwrite_dfns = 100 diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml new file mode 100644 index 000000000..9f0a16f23 --- /dev/null +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml @@ -0,0 +1,71 @@ +# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions. +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +T_e = 1.0 +T_wall = 1.0 +initial_density1 = 0.5 +initial_temperature1 = 1.0 +initial_density2 = 0.5 +initial_temperature2 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.001 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +z_IC_option2 = "sinusoid" +z_IC_density_amplitude2 = 0.001 +z_IC_density_phase2 = 0.0 +z_IC_upar_amplitude2 = 0.0 +z_IC_upar_phase2 = 0.0 +z_IC_temperature_amplitude2 = 0.0 +z_IC_temperature_phase2 = 0.0 +vpa_IC_option1 = "isotropic-beam" +#vpa_IC_option1 = "directed-beam" +#vpa_IC_v01 = +#vpa_IC_vth0 = +#vpa_IC_vpa0 = +#vpa_IC_vperp0 = +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false +z_ngrid = 1 +z_nelement = 1 +z_nelement_local = 1 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 6.0 +vpa_bc = "zero" +vpa_discretization = "gausslegendre_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 3.0 +vperp_discretization = "gausslegendre_pseudospectral" +# Fokker-Planck operator requires the "gausslegendre_pseudospectral +# options for the vpa and vperp grids + +[fokker_planck_collisions] +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.0 +slowing_down_test = true +frequency_option = "manual" +use_fokker_planck = true + +[timestepping] +nstep = 1000 +dt = 1.0e-5 +nwrite = 50 +nwrite_dfns = 50 diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml new file mode 100644 index 000000000..9fe59dbe2 --- /dev/null +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml @@ -0,0 +1,65 @@ +# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions. +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +T_e = 1.0 +T_wall = 1.0 +initial_density1 = 0.5 +initial_temperature1 = 1.0 +initial_density2 = 0.5 +initial_temperature2 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.001 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +z_IC_option2 = "sinusoid" +z_IC_density_amplitude2 = 0.001 +z_IC_density_phase2 = 0.0 +z_IC_upar_amplitude2 = 0.0 +z_IC_upar_phase2 = 0.0 +z_IC_temperature_amplitude2 = 0.0 +z_IC_temperature_phase2 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false +z_ngrid = 1 +z_nelement = 1 +z_nelement_local = 1 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 6.0 +vpa_bc = "zero" +vpa_discretization = "gausslegendre_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 3.0 +vperp_discretization = "gausslegendre_pseudospectral" +# Fokker-Planck operator requires the "gausslegendre_pseudospectral +# options for the vpa and vperp grids + +[fokker_planck_collisions] +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.0 +slowing_down_test = true +frequency_option = "manual" +use_fokker_planck = true + +[timestepping] +nstep = 1000 +dt = 1.0e-5 +nwrite = 50 +nwrite_dfns = 50 diff --git a/examples/fokker-planck/fokker-planck-relaxation.toml b/examples/fokker-planck/fokker-planck-relaxation.toml index cea04911c..f2d5651a1 100644 --- a/examples/fokker-planck/fokker-planck-relaxation.toml +++ b/examples/fokker-planck/fokker-planck-relaxation.toml @@ -63,4 +63,3 @@ nstep = 5000 dt = 1.0e-2 nwrite = 5000 nwrite_dfns = 5000 -split_operators = false diff --git a/moment_kinetics/src/communication.jl b/moment_kinetics/src/communication.jl index 2550d0358..ee1b84f3a 100644 --- a/moment_kinetics/src/communication.jl +++ b/moment_kinetics/src/communication.jl @@ -120,7 +120,7 @@ function __init__() if !MPI.Initialized() MPI.Init() end - + comm_world.val = MPI.COMM_WORLD.val global_rank[] = MPI.Comm_rank(comm_world) diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl index dae38e039..95ce99da6 100644 --- a/moment_kinetics/src/coordinates.jl +++ b/moment_kinetics/src/coordinates.jl @@ -192,7 +192,8 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing cell_width, igrid, ielement, imin, imax, igrid_full, input.discretization, input.fd_option, input.cheb_option, input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch), scratch_shared, scratch_shared2, scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm, - local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option, element_boundaries) + local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option, + element_boundaries) if coord.n == 1 && occursin("v", coord.name) spectral = null_velocity_dimension_info() diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 012183483..0352f9a1f 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -33,6 +33,7 @@ export init_fokker_planck_collisions, fokkerplanck_arrays_struct export init_fokker_planck_collisions_weak_form export explicit_fokker_planck_collisions_weak_form! export explicit_fokker_planck_collisions! +export explicit_fp_collisions_weak_form_Maxwellian_cross_species! export calculate_Maxwellian_Rosenbluth_coefficients export get_local_Cssp_coefficients!, init_fokker_planck_collisions # testing @@ -87,8 +88,15 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) # begin default inputs (as kwargs) use_fokker_planck = false, nuii = -1.0, - frequency_option = "reference_parameters") - + frequency_option = "reference_parameters", + self_collisions = true, + slowing_down_test = false, + sd_density = 1.0, + sd_temp = 0.01, + sd_q = 1.0, + sd_mi = 0.25, + sd_me = 0.25/1836.0, + Zi = 1.0) # ensure that the collision frequency is consistent with the input option frequency_option = input_section["frequency_option"] if frequency_option == "reference_parameters" @@ -107,6 +115,30 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) end input = Dict(Symbol(k)=>v for (k,v) in input_section) #println(input) + if input_section["slowing_down_test"] + # calculate nu_alphae and vc3 (critical speed of slowing down) + # as a diagnostic to aid timestep choice + # nu_alphae/nuref + Zalpha = input_section["Zi"] + Zi = input_section["sd_q"] + ni = input_section["sd_density"] + ne = Zi*ni+Zalpha # assume unit density for alphas in initial state + Te = input_section["sd_temp"] + me = input_section["sd_me"] + mi = input_section["sd_mi"] + nu_alphae = (8.0/(3.0*sqrt(pi)))*ne*sqrt(me)*(Te^(-1.5))*(Zalpha^(2.0)) + vc3 = (3.0/4.0)*sqrt(pi)*(Zi^2)*(ni/ne)*(1.0/mi)*(1.0/sqrt(me))*(Te^(1.5)) + nu_alphaalpha = Zalpha^4 + nu_alphai = nu_alphae*vc3 + if global_rank[] == 0 + println("slowing_down_test = true") + println("nu_alphaalpha/nuref = $nu_alphaalpha") + println("nu_alphai/nuref = $nu_alphai") + println("nu_alphae/nuref = $nu_alphae") + println("vc3/cref^3 = $vc3") + println("critical speed vc/cref = ",vc3^(1.0/3.0)) + end + end return fkpl_collisions_input(; input...) end @@ -114,7 +146,7 @@ end # begin functions associated with the weak-form operator # where the potentials are computed by elliptic solve ######################################################## - + """ function that initialises the arrays needed for Fokker Planck collisions using numerical integration to compute the Rosenbluth potentials only @@ -199,6 +231,71 @@ function init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_sp return fka end +""" +Function for advancing with the explicit, weak-form, self-collision operator +using the existing method for computing the Rosenbluth potentials, with +the addition of cross-species collisions against fixed Maxwellian distribution functions +""" +function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_in,dSdt,composition,collisions,dt, + fkpl_arrays::fokkerplanck_weakform_arrays_struct, + r, z, vperp, vpa, vperp_spectral, vpa_spectral; + diagnose_entropy_production=false) + # N.B. only self-collisions are currently supported + # This can be modified by adding a loop over s' below + n_ion_species = composition.n_ion_species + @boundscheck vpa.n == size(pdf_out,1) || throw(BoundsError(pdf_out)) + @boundscheck vperp.n == size(pdf_out,2) || throw(BoundsError(pdf_out)) + @boundscheck z.n == size(pdf_out,3) || throw(BoundsError(pdf_out)) + @boundscheck r.n == size(pdf_out,4) || throw(BoundsError(pdf_out)) + @boundscheck n_ion_species == size(pdf_out,5) || throw(BoundsError(pdf_out)) + @boundscheck vpa.n == size(pdf_in,1) || throw(BoundsError(pdf_in)) + @boundscheck vperp.n == size(pdf_in,2) || throw(BoundsError(pdf_in)) + @boundscheck z.n == size(pdf_in,3) || throw(BoundsError(pdf_in)) + @boundscheck r.n == size(pdf_in,4) || throw(BoundsError(pdf_in)) + @boundscheck n_ion_species == size(pdf_in,5) || throw(BoundsError(pdf_in)) + @boundscheck z.n == size(dSdt,1) || throw(BoundsError(dSdt)) + @boundscheck r.n == size(dSdt,2) || throw(BoundsError(dSdt)) + @boundscheck n_ion_species == size(dSdt,3) || throw(BoundsError(dSdt)) + + fkin = collisions.fkpl + # masses charge numbers and collision frequencies + mref = 1.0 # generalise if multiple ions evolved + Zs = fkin.Zi # generalise if multiple ions evolved + nuref = fkin.nuii # generalise if multiple ions evolved + msp = [fkin.sd_mi, fkin.sd_me] + Zsp = [fkin.sd_q, -1.0] + # assume here that ne = sum_i n_i and that initial condition + # for beam ions has unit density + ns = 1.0 # initial density of evolved ions (hardcode to be the same as initial conditions) + # get electron density from quasineutrality ne = sum_s Zs ns + densp = [fkin.sd_density, fkin.sd_q*fkin.sd_density+ns*Zs] + uparsp = [0.0, 0.0] + vthsp = [sqrt(fkin.sd_temp/msp[1]), sqrt(fkin.sd_temp/msp[2])] + + # N.B. parallelisation using special 'anyv' region + begin_s_r_z_anyv_region() + @loop_s_r_z is ir iz begin + # computes sum over s' of C[Fs,Fs'] with Fs' an assumed Maxwellian + @views fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!(pdf_in[:,:,iz,ir,is], + nuref,mref,Zs,msp,Zsp,densp,uparsp,vthsp, + fkpl_arrays,vperp,vpa,vperp_spectral,vpa_spectral) + # 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 sure that the cross-species terms conserve density + density_conserving_correction!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, + fkpl_arrays.S_dummy) + # advance this part of s,r,z with the resulting sum_s' C[Fs,Fs'] + begin_anyv_vperp_vpa_region() + CC = fkpl_arrays.CC + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] += dt*CC[ivpa,ivperp] + end + + end + return nothing +end + + """ Function for advancing with the explicit, weak-form, self-collision operator """ @@ -226,7 +323,9 @@ function explicit_fokker_planck_collisions_weak_form!(pdf_out,pdf_in,dSdt,compos # masses and collision frequencies ms, msp = 1.0, 1.0 # generalise! - nussp = collisions.fkpl.nuii # generalise! + nuref = collisions.fkpl.nuii # generalise! + Zi = collisions.fkpl.Zi # generalise! + nussp = nuref*(Zi^4) # include charge number factor for self collisions # N.B. parallelisation using special 'anyv' region begin_s_r_z_anyv_region() @loop_s_r_z is ir iz begin @@ -267,12 +366,18 @@ Function for evaluating \$C_{ss'} = C_{ss'}[F_s,F_{s'}]\$ The result is stored in the array `fkpl_arrays.CC`. -The normalised collision frequency is defined by +The normalised collision frequency for collisions between species s and s' is defined by ```math \\tilde{\\nu}_{ss'} = \\frac{L_{\\mathrm{ref}}}{c_{\\mathrm{ref}}}\\frac{\\gamma_{ss'} n_\\mathrm{ref}}{m_s^2 c_\\mathrm{ref}^3} ``` with \$\\gamma_{ss'} = 2 \\pi (Z_s Z_{s'})^2 e^4 \\ln \\Lambda_{ss'} / (4 \\pi \\epsilon_0)^2\$. +The input parameter to this code is +```math +\\tilde{\\nu}_{ii} = \\frac{L_{\\mathrm{ref}}}{c_{\\mathrm{ref}}}\\frac{\\gamma_\\mathrm{ref} n_\\mathrm{ref}}{m_\\mathrm{ref}^2 c_\\mathrm{ref}^3} +``` +with \$\\gamma_\\mathrm{ref} = 2 \\pi e^4 \\ln \\Lambda_{ii} / (4 \\pi +\\epsilon_0)^2\$. This means that \$\\tilde{\\nu}_{ss'} = (Z_s Z_{s'})^2\\tilde{\\nu}_\\mathrm{ref}\$ and this conversion is handled explicitly in the code with the charge number input provided by the user. """ function fokker_planck_collision_operator_weak_form!(ffs_in,ffsp_in,ms,msp,nussp, fkpl_arrays::fokkerplanck_weakform_arrays_struct, @@ -383,6 +488,98 @@ function fokker_planck_collision_operator_weak_form!(ffs_in,ffsp_in,ms,msp,nussp return nothing end +""" +Function for computing the collision operator +```math +\\sum_{s^\\prime} C[F_{s},F_{s^\\prime}] +``` +when +```math +F_{s^\\prime} +``` +is an analytically specified Maxwellian distribution +""" +function fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!(ffs_in, + nuref::mk_float,ms::mk_float,Zs::mk_float, + msp::Array{mk_float,1},Zsp::Array{mk_float,1}, + densp::Array{mk_float,1}, + uparsp::Array{mk_float,1},vthsp::Array{mk_float,1}, + fkpl_arrays::fokkerplanck_weakform_arrays_struct, + vperp, vpa, vperp_spectral, vpa_spectral) + @boundscheck vpa.n == size(ffs_in,1) || throw(BoundsError(ffs_in)) + @boundscheck vperp.n == size(ffs_in,2) || throw(BoundsError(ffs_in)) + + # extract the necessary precalculated and buffer arrays from fokkerplanck_arrays + rhsvpavperp = fkpl_arrays.rhsvpavperp + lu_obj_MM = fkpl_arrays.lu_obj_MM + YY_arrays = fkpl_arrays.YY_arrays + + CC = fkpl_arrays.CC + GG = fkpl_arrays.GG + HH = fkpl_arrays.HH + dHdvpa = fkpl_arrays.dHdvpa + dHdvperp = fkpl_arrays.dHdvperp + dGdvperp = fkpl_arrays.dGdvperp + d2Gdvperp2 = fkpl_arrays.d2Gdvperp2 + d2Gdvpa2 = fkpl_arrays.d2Gdvpa2 + d2Gdvperpdvpa = fkpl_arrays.d2Gdvperpdvpa + FF = fkpl_arrays.FF + dFdvpa = fkpl_arrays.dFdvpa + dFdvperp = fkpl_arrays.dFdvperp + + # number of primed species + nsp = size(msp,1) + + begin_anyv_vperp_vpa_region() + # fist set dummy arrays for coefficients to zero + @loop_vperp_vpa ivperp ivpa begin + d2Gdvpa2[ivpa,ivperp] = 0.0 + d2Gdvperp2[ivpa,ivperp] = 0.0 + d2Gdvperpdvpa[ivpa,ivperp] = 0.0 + dHdvpa[ivpa,ivperp] = 0.0 + dHdvperp[ivpa,ivperp] = 0.0 + end + # sum the contributions from the potentials, including order unity factors that differ between species + # making use of the Linearity of the operator in Fsp + # note that here we absorb ms/msp and Zsp^2 into the definition of the potentials, and we pass + # ms = msp = 1 to the collision operator assembly routine so that we can use a single array to include + # the contribution to the summed Rosenbluth potential from all the species + for isp in 1:nsp + dens = densp[isp] + upar = uparsp[isp] + vth = vthsp[isp] + ZZ = (Zsp[isp]*Zs)^2 # factor from gamma_ss' + @loop_vperp_vpa ivperp ivpa begin + d2Gdvpa2[ivpa,ivperp] += ZZ*d2Gdvpa2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + d2Gdvperp2[ivpa,ivperp] += ZZ*d2Gdvperp2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + d2Gdvperpdvpa[ivpa,ivperp] += ZZ*d2Gdvperpdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + dHdvpa[ivpa,ivperp] += ZZ*(ms/msp[isp])*dHdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + dHdvperp[ivpa,ivperp] += ZZ*(ms/msp[isp])*dHdvperp_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + end + end + # Need to synchronize as these arrays may be read outside the locally-owned set of + # ivperp, ivpa indices in assemble_explicit_collision_operator_rhs_parallel!() + # assemble the RHS of the collision operator matrix eq + + _anyv_subblock_synchronize() + assemble_explicit_collision_operator_rhs_parallel!(rhsvpavperp,ffs_in, + d2Gdvpa2,d2Gdvperpdvpa,d2Gdvperp2, + dHdvpa,dHdvperp,1.0,1.0,nuref, + vpa,vperp,YY_arrays) + + # solve the collision operator matrix eq + begin_anyv_region() + @anyv_serial_region begin + # sc and rhsc are 1D views of the data in CC and rhsc, created so that we can use + # the 'matrix solve' functionality of ldiv!() from the LinearAlgebra package + sc = vec(CC) + rhsc = vec(rhsvpavperp) + # invert mass matrix and fill fc + ldiv!(sc, lu_obj_MM, rhsc) + end + return nothing +end + # solves A x = b for a matrix of the form # A00 0 A02 # 0 A11 A12 @@ -478,6 +675,36 @@ function conserving_corrections!(CC,pdf_in,vpa,vperp,dummy_vpavperp) end end +function density_conserving_correction!(CC,pdf_in,vpa,vperp,dummy_vpavperp) + begin_anyv_region() + x0 = 0.0 + @anyv_serial_region begin + # In principle the integrations here could be shared among the processes in the + # 'anyv' subblock, but this block is not a significant part of the cost of the + # collision operator, so probably not worth the complication. + + # compute density of the input pdf + dens = get_density(pdf_in, vpa, vperp) + + # compute density of the numerical collision operator + dn = get_density(CC, vpa, vperp) + + # obtain the coefficient for the correction + x0 = dn/dens + end + + # Broadcast x0 to all processes in the 'anyv' subblock + param_vec = [x0] + MPI.Bcast!(param_vec, 0, comm_anyv_subblock[]) + x0 = param_vec[1] + + # correct CC + begin_anyv_vperp_vpa_region() + @loop_vperp_vpa ivperp ivpa begin + CC[ivpa,ivperp] -= x0*pdf_in[ivpa,ivperp] + end +end + ###################################################### # end functions associated with the weak-form operator diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index 9c3e9a03d..ca8697628 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -29,6 +29,7 @@ export enforce_zero_bc! export allocate_rosenbluth_potential_boundary_data export calculate_rosenbluth_potential_boundary_data_exact! export test_rosenbluth_potential_boundary_data +export interpolate_2D_vspace! # Import moment_kinetics so that we can refer to it in docstrings import moment_kinetics @@ -2317,4 +2318,97 @@ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral) end end +""" +function to interpolate f(vpa,vperp) from one +velocity grid to another, assuming that both +grids are represented by vpa, vperp in normalised units, +but have different normalisation factors +defining the meaning of these grids in physical units. + +E.g. vpai, vperpi = ci * vpa, ci * vperp + vpae, vperpe = ce * vpa, ce * vperp + +with ci = sqrt(Ti/mi), ce = sqrt(Te/mi) + +scalefac = ci / ce is the ratio of the +two reference speeds + +""" +function interpolate_2D_vspace!(pdf_out,pdf_in,vpa,vperp,scalefac) + + begin_anyv_vperp_vpa_region() + # loop over points in the output interpolated dataset + @loop_vperp ivperp begin + vperp_val = vperp.grid[ivperp]*scalefac + # get element for interpolation data + iel_vperp = ielement_loopup(vperp_val,vperp) + if iel_vperp < 1 # vperp_interp outside of range of vperp.grid + @loop_vpa ivpa begin + pdf_out[ivpa,ivperp] = 0.0 + end + continue + else + # get nodes for interpolation + ivperpmin, ivperpmax = vperp.igrid_full[1,iel_vperp], vperp.igrid_full[vperp.ngrid,iel_vperp] + vperp_nodes = vperp.grid[ivperpmin:ivperpmax] + #print("vperp: ",iel_vperp, " ", vperp_nodes," ",vperp_val) + + end + @loop_vpa ivpa begin + vpa_val = vpa.grid[ivpa]*scalefac + # get element for interpolation data + iel_vpa = ielement_loopup(vpa_val,vpa) + if iel_vpa < 1 # vpa_interp outside of range of vpa.grid + pdf_out[ivpa,ivperp] = 0.0 + continue + else + # get nodes for interpolation + ivpamin, ivpamax = vpa.igrid_full[1,iel_vpa], vpa.igrid_full[vpa.ngrid,iel_vpa] + vpa_nodes = vpa.grid[ivpamin:ivpamax] + #print("vpa: ", iel_vpa, " ", vpa_nodes," ",vpa_val) + + # do the interpolation + pdf_out[ivpa,ivperp] = 0.0 + for ivperpgrid in 1:vperp.ngrid + # index for referencing pdf_in on orginal grid + ivperpp = vperp.igrid_full[ivperpgrid,iel_vperp] + # interpolating polynomial value at ivperpp for interpolation + vperppoly = lagrange_poly(ivperpgrid,vperp_nodes,vperp_val) + for ivpagrid in 1:vpa.ngrid + # index for referencing pdf_in on orginal grid + ivpap = vpa.igrid_full[ivpagrid,iel_vpa] + # interpolating polynomial value at ivpap for interpolation + vpapoly = lagrange_poly(ivpagrid,vpa_nodes,vpa_val) + pdf_out[ivpa,ivperp] += vpapoly*vperppoly*pdf_in[ivpap,ivperpp] + end + end + end + end + end + return nothing +end + +""" +function to find the element in which x sits +""" +function ielement_loopup(x,coord) + xebs = coord.element_boundaries + nelement = coord.nelement_global + zero = 1.0e-14 + ielement = -1 + # find the element + for j in 1:nelement + # check for internal points + if (x - xebs[j])*(xebs[j+1] - x) > zero + ielement = j + break + # check for boundary points + elseif (abs(x-xebs[j]) < 100*zero) || (abs(x-xebs[j+1]) < 100*zero && j == nelement) + ielement = j + break + end + end + return ielement +end + end diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index bc82f058d..7f724c6bc 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -665,6 +665,31 @@ function init_ion_pdf_over_density!(pdf, spec, composition, vpa, vperp, z, @loop_z_vperp iz ivperp begin @. pdf[:,ivperp,iz] = (vpa.grid + 0.5*vpa.L)^spec.vpa_IC.monomial_degree end + elseif spec.vpa_IC.initialization_option == "isotropic-beam" + v0 = spec.vpa_IC.v0 #0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2) # birth speed of beam + vth0 = spec.vpa_IC.vth0 + v4norm = (v0^2)*(vth0^2) # spread of the beam in speed is vth0 + @loop_z iz begin + @loop_vperp_vpa ivperp ivpa begin + v2 = (vpa.grid[ivpa])^2 + vperp.grid[ivperp]^2 - v0^2 + pdf[ivpa,ivperp,iz] = exp(-(v2^2)/v4norm) + end + normfac = integrate_over_vspace(view(pdf,:,:,iz), vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts) + @. pdf[:,:,iz] /= normfac + end + elseif spec.vpa_IC.initialization_option == "directed-beam" + vpa0 = spec.vpa_IC.vpa0 #0.25*0.5*abs(vpa.L) # centre of beam in vpa + vperp0 = spec.vpa_IC.vperp0 #0.5*abs(vperp.L) # centre of beam in vperp + vth0 = spec.vpa_IC.vth0 #0.05*sqrt(vperp.L^2 + (0.5*vpa.L)^2) # width of beam in v + @loop_z iz begin + @loop_vperp_vpa ivperp ivpa begin + v2 = (vpa.grid[ivpa] - vpa0)^2 + (vperp.grid[ivperp] - vperp0)^2 + v2norm = vth0^2 + pdf[ivpa,ivperp,iz] = exp(-v2/v2norm) + end + normfac = integrate_over_vspace(view(pdf,:,:,iz), vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts) + @. pdf[:,:,iz] /= normfac + end end return nothing end diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 690776b67..fe7661e9a 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -239,6 +239,11 @@ struct initial_condition_input temperature_phase::mk_float # inputs for "monomial" initial condition monomial_degree::mk_int + # inputs for "isotropic-beam", "directed-beam" initial conditions + v0::mk_float + vth0::mk_float + vpa0::mk_float + vperp0::mk_float end """ @@ -351,13 +356,34 @@ Base.@kwdef struct krook_collisions_input end Base.@kwdef struct fkpl_collisions_input + # option to check if fokker planck frequency should be > 0 use_fokker_planck::Bool - # ion-ion self collision frequency - # nu_{ss'} = gamma_{ss'} n_{ref} / (m_s)^2 (c_{ref})^3 - # with gamma_ss' = 2 pi (Z_s Z_s')^2 e^4 ln \Lambda_{ss'} / (4 pi \epsilon_0)^2 + # ion-ion self collision frequency (for a species with Z = 1) + # nu_{ii} = (L/c_{ref}) * gamma_{ref} n_{ref} /(m_s)^2 (c_{ref})^3 + # with gamma_ref = 2 pi e^4 ln \Lambda_{ii} / (4 pi \epsilon_0)^2 + # and ln \Lambda_{ii} the Coulomb logarithm for ion-ion collisions nuii::mk_float + # option to determine if self collisions are used (for physics test) + self_collisions::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 frequency_option::String # "manual" # "reference_parameters" + # options for fixed Maxwellian species in slowing down test operator + # ion density - electron density determined from quasineutrality + sd_density::mk_float + # ion temperature - electron temperature assumed identical + sd_temp::mk_float + # ion charge number of fixed Maxwellian species + sd_q::mk_float + # ion mass with respect to reference + sd_mi::mk_float + # electron mass with respect to reference + sd_me::mk_float + # charge number of evolved ion species + # kept here because charge number different from 1 + # is not supported for other physics features + Zi::mk_float end """ diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 456779ec7..56358a27b 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -535,19 +535,23 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) species.ion[is].z_IC.density_amplitude, species.ion[is].z_IC.density_phase, species.ion[is].z_IC.upar_amplitude, species.ion[is].z_IC.upar_phase, species.ion[is].z_IC.temperature_amplitude, species.ion[is].z_IC.temperature_phase, - species.ion[is].z_IC.monomial_degree) + species.ion[is].z_IC.monomial_degree, 0.0, 0.0, 0.0, 0.0) r_IC = initial_condition_input(species.ion[is].r_IC.initialization_option, species.ion[is].r_IC.width, species.ion[is].r_IC.wavenumber, species.ion[is].r_IC.density_amplitude, species.ion[is].r_IC.density_phase, species.ion[is].r_IC.upar_amplitude, species.ion[is].r_IC.upar_phase, species.ion[is].r_IC.temperature_amplitude, species.ion[is].r_IC.temperature_phase, - species.ion[is].r_IC.monomial_degree) + species.ion[is].r_IC.monomial_degree, 0.0, 0.0, 0.0, 0.0) vpa_IC = initial_condition_input(species.ion[is].vpa_IC.initialization_option, species.ion[is].vpa_IC.width, species.ion[is].vpa_IC.wavenumber, species.ion[is].vpa_IC.density_amplitude, species.ion[is].vpa_IC.density_phase, species.ion[is].vpa_IC.upar_amplitude, species.ion[is].vpa_IC.upar_phase, species.ion[is].vpa_IC.temperature_amplitude, - species.ion[is].vpa_IC.temperature_phase, species.ion[is].vpa_IC.monomial_degree) + species.ion[is].vpa_IC.temperature_phase, species.ion[is].vpa_IC.monomial_degree, + get(scan_input, "vpa_IC_v0_$is", 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), + get(scan_input, "vpa_IC_vth0$is", 0.1*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), + get(scan_input, "vpa_IC_vpa0$is", 0.25*0.5*abs(vpa.L)), + get(scan_input, "vpa_IC_vperp0$is", 0.5*abs(vperp.L))) species_ion_immutable[is] = species_parameters(species_type, species.ion[is].initial_temperature, species.ion[is].initial_density, z_IC, r_IC, vpa_IC) end @@ -559,19 +563,23 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) species.neutral[is].z_IC.density_amplitude, species.neutral[is].z_IC.density_phase, species.neutral[is].z_IC.upar_amplitude, species.neutral[is].z_IC.upar_phase, species.neutral[is].z_IC.temperature_amplitude, species.neutral[is].z_IC.temperature_phase, - species.neutral[is].z_IC.monomial_degree) + species.neutral[is].z_IC.monomial_degree, 0.0, 0.0, 0.0, 0.0) r_IC = initial_condition_input(species.neutral[is].r_IC.initialization_option, species.neutral[is].r_IC.width, species.neutral[is].r_IC.wavenumber, species.neutral[is].r_IC.density_amplitude, species.neutral[is].r_IC.density_phase, species.neutral[is].r_IC.upar_amplitude, species.neutral[is].r_IC.upar_phase, species.neutral[is].r_IC.temperature_amplitude, species.neutral[is].r_IC.temperature_phase, - species.neutral[is].r_IC.monomial_degree) + species.neutral[is].r_IC.monomial_degree, 0.0, 0.0, 0.0, 0.0) vpa_IC = initial_condition_input(species.neutral[is].vpa_IC.initialization_option, species.neutral[is].vpa_IC.width, species.neutral[is].vpa_IC.wavenumber, species.neutral[is].vpa_IC.density_amplitude, species.neutral[is].vpa_IC.density_phase, species.neutral[is].vpa_IC.upar_amplitude, species.neutral[is].vpa_IC.upar_phase, species.neutral[is].vpa_IC.temperature_amplitude, - species.neutral[is].vpa_IC.temperature_phase, species.neutral[is].vpa_IC.monomial_degree) + species.neutral[is].vpa_IC.temperature_phase, species.neutral[is].vpa_IC.monomial_degree, + get(scan_input, "vpa_IC_v0$is", 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), + get(scan_input, "vpa_IC_vth0$is", 0.1*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), + get(scan_input, "vpa_IC_vpa0$is", 0.25*0.5*abs(vpa.L)), + get(scan_input, "vpa_IC_vperp0$is", 0.5*abs(vperp.L))) species_neutral_immutable[is] = species_parameters(species_type, species.neutral[is].initial_temperature, species.neutral[is].initial_density, z_IC, r_IC, vpa_IC) end diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 157257955..e5cdc0ef4 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -54,6 +54,7 @@ using ..force_balance: force_balance!, neutral_force_balance! using ..energy_equation: energy_equation!, neutral_energy_equation! using ..em_fields: setup_em_fields, update_phi! using ..fokker_planck: init_fokker_planck_collisions_weak_form, explicit_fokker_planck_collisions_weak_form! +using ..fokker_planck: explicit_fp_collisions_weak_form_Maxwellian_cross_species! using ..gyroaverages: init_gyro_operators, gyroaverage_pdf! using ..manufactured_solns: manufactured_sources using ..advection: advection_info @@ -2191,9 +2192,18 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, # advance with the Fokker-Planck self-collision operator if advance.explicit_weakform_fp_collisions update_entropy_diagnostic = (istage == 1) - explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt,composition,collisions,dt, - fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, + if collisions.fkpl.self_collisions + # self collisions for each species + explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt,composition, + collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, + diagnose_entropy_production = update_entropy_diagnostic) + end + if collisions.fkpl.slowing_down_test + # include cross-collsions with fixed Maxwellian backgrounds + explicit_fp_collisions_weak_form_Maxwellian_cross_species!(fvec_out.pdf,fvec_in.pdf,moments.ion.dSdt, + composition,collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral, diagnose_entropy_production = update_entropy_diagnostic) + end end # End of advance for distribution function diff --git a/moment_kinetics/test/fokker_planck_tests.jl b/moment_kinetics/test/fokker_planck_tests.jl index 54266f518..edfd0d4d4 100644 --- a/moment_kinetics/test/fokker_planck_tests.jl +++ b/moment_kinetics/test/fokker_planck_tests.jl @@ -15,6 +15,7 @@ using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_ppe 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 +using moment_kinetics.fokker_planck: density_conserving_correction!, fokker_planck_collision_operator_weak_form_Maxwellian_Fsp! using moment_kinetics.fokker_planck_test: print_test_data, fkpl_error_data, allocate_error_data #, plot_test_data using moment_kinetics.fokker_planck_test: F_Maxwellian, G_Maxwellian, H_Maxwellian using moment_kinetics.fokker_planck_test: d2Gdvpa2_Maxwellian, d2Gdvperp2_Maxwellian, d2Gdvperpdvpa_Maxwellian, dGdvperp_Maxwellian @@ -22,6 +23,7 @@ using moment_kinetics.fokker_planck_test: dHdvperp_Maxwellian, dHdvpa_Maxwellian using moment_kinetics.fokker_planck_calculus: calculate_rosenbluth_potentials_via_elliptic_solve!, calculate_rosenbluth_potential_boundary_data_exact! using moment_kinetics.fokker_planck_calculus: test_rosenbluth_potential_boundary_data, allocate_rosenbluth_potential_boundary_data using moment_kinetics.fokker_planck_calculus: enforce_vpavperp_BCs!, calculate_rosenbluth_potentials_via_direct_integration! +using moment_kinetics.fokker_planck_calculus: interpolate_2D_vspace! function create_grids(ngrid,nelement_vpa,nelement_vperp; Lvpa=12.0,Lvperp=6.0) @@ -70,6 +72,80 @@ function runtests() @testset "Fokker Planck tests" verbose=use_verbose begin println("Fokker Planck tests") + @testset " - test Lagrange-polynomial 2D interpolation" begin + println(" - test Lagrange-polynomial 2D interpolation") + ngrid = 9 + nelement_vpa = 16 + nelement_vperp = 8 + vpa, vpa_spectral, vperp, vperp_spectral = create_grids(ngrid,nelement_vpa,nelement_vperp, + Lvpa=8.0,Lvperp=4.0) + + # electron pdf on electron grids + Fe = allocate_shared_float(vpa.n,vperp.n) + # electron pdf on ion normalised grids + Fe_interp_ion_units = allocate_shared_float(vpa.n,vperp.n) + # exact value for comparison + Fe_exact_ion_units = allocate_shared_float(vpa.n,vperp.n) + # ion pdf on ion grids + Fi = allocate_shared_float(vpa.n,vperp.n) + # ion pdf on electron normalised grids + Fi_interp_electron_units = allocate_shared_float(vpa.n,vperp.n) + # exact value for comparison + Fi_exact_electron_units = allocate_shared_float(vpa.n,vperp.n) + # test array + F_err = allocate_float(vpa.n,vperp.n) + + dense = 1.0 + upare = 0.0 # upare in electron reference units + vthe = 1.0 # vthe in electron reference units + densi = 1.0 + upari = 0.0 # upari in ion reference units + vthi = 1.0 # vthi in ion reference units + # reference speeds for electrons and ions + cref_electron = 60.0 + cref_ion = 1.0 + # scale factor for change of reference speed + scalefac = cref_ion/cref_electron + + begin_serial_region() + @serial_region begin + @loop_vperp_vpa ivperp ivpa begin + Fe[ivpa,ivperp] = F_Maxwellian(dense,upare,vthe,vpa,vperp,ivpa,ivperp) + Fe_exact_ion_units[ivpa,ivperp] = F_Maxwellian(dense,upare/scalefac,vthe/scalefac,vpa,vperp,ivpa,ivperp)/(scalefac^3) + Fi[ivpa,ivperp] = F_Maxwellian(densi,upari,vthi,vpa,vperp,ivpa,ivperp) + Fi_exact_electron_units[ivpa,ivperp] = (scalefac^3)*F_Maxwellian(densi,upari*scalefac,vthi*scalefac,vpa,vperp,ivpa,ivperp) + end + end + + begin_s_r_z_anyv_region() + interpolate_2D_vspace!(Fe_interp_ion_units,Fe,vpa,vperp,scalefac) + #println("Fe",Fe) + #println("Fe interp",Fe_interp_ion_units) + #println("Fe exact",Fe_exact_ion_units) + interpolate_2D_vspace!(Fi_interp_electron_units,Fi,vpa,vperp,1.0/scalefac) + #println("Fi",Fi) + #println("Fi interp", Fi_interp_electron_units) + #println("Fi exact",Fi_exact_electron_units) + + begin_serial_region() + # check the result + @serial_region begin + # for electron data on ion grids + @. F_err = abs(Fe_interp_ion_units - Fe_exact_ion_units) + max_F_err = maximum(F_err) + max_F = maximum(Fe_exact_ion_units) + #println(max_F) + @test max_F_err < 3.0e-8 * max_F + # for ion data on electron grids + @. F_err = abs(Fi_interp_electron_units - Fi_exact_electron_units) + max_F_err = maximum(F_err) + max_F = maximum(Fi_exact_electron_units) + #println(max_F) + @test max_F_err < 3.0e-8 * max_F + end + + end + @testset " - test weak-form 2D differentiation" begin # tests the correct definition of mass and stiffness matrices in 2D println(" - test weak-form 2D differentiation") @@ -142,7 +218,6 @@ function runtests() nelement_vperp = 4 vpa, vpa_spectral, vperp, vperp_spectral = create_grids(ngrid,nelement_vpa,nelement_vperp, Lvpa=12.0,Lvperp=6.0) - nc_global = vpa.n*vperp.n begin_serial_region() fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral, precompute_weights=true, print_to_screen=print_to_screen) @@ -290,7 +365,6 @@ function runtests() nelement_vperp = 4 vpa, vpa_spectral, vperp, vperp_spectral = create_grids(ngrid,nelement_vpa,nelement_vperp, Lvpa=12.0,Lvperp=6.0) - nc_global = vpa.n*vperp.n begin_serial_region() fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral, precompute_weights=true, print_to_screen=print_to_screen) @@ -440,6 +514,102 @@ function runtests() finalize_comms!() end + @testset " - test weak-form (slowing-down) collision operator calculation" begin + println(" - test weak-form (slowing-down) collision operator calculation") + ngrid = 9 + nelement_vpa = 16 + nelement_vperp = 8 + vpa, vpa_spectral, vperp, vperp_spectral = create_grids(ngrid,nelement_vpa,nelement_vperp, + Lvpa=12.0,Lvperp=6.0) + begin_serial_region() + fkpl_arrays = init_fokker_planck_collisions_weak_form(vpa,vperp,vpa_spectral,vperp_spectral, + precompute_weights=true, print_to_screen=print_to_screen) + + @testset "slowing_down_test=true test_numerical_conserving_terms=$test_numerical_conserving_terms" for test_numerical_conserving_terms in (true,false) + + dummy_array = allocate_float(vpa.n,vperp.n) + Fs_M = allocate_float(vpa.n,vperp.n) + F_M = allocate_float(vpa.n,vperp.n) + C_M_num = allocate_shared_float(vpa.n,vperp.n) + C_M_exact = allocate_float(vpa.n,vperp.n) + C_M_err = allocate_float(vpa.n,vperp.n) + + # pick a set of parameters that represent slowing down + # on slow ions and faster electrons, but which are close + # enough to 1 for errors comparable to the self-collision operator + # increasing or reducing vth, mass increases the errors + dens, upar, vth = 1.0, 1.0, 1.0 + mref = 1.0 + Zref = 1.0 + msp = [1.0,0.2]#[0.25, 0.25/1836.0] + Zsp = [0.5,0.5]#[0.5, 0.5] + denssp = [1.0,1.0]#[1.0, 1.0] + uparsp = [0.0,0.0]#[0.0, 0.0] + vthsp = [sqrt(0.5/msp[1]), sqrt(0.5/msp[2])]#[sqrt(0.01/msp[1]), sqrt(0.01/msp[2])] + nsprime = size(msp,1) + nuref = 1.0 + + begin_serial_region() + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + Fs_M[ivpa,ivperp] = F_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + C_M_exact[ivpa,ivperp] = 0.0 + end + end + # sum up contributions to cross-collision operator + for isp in 1:nsprime + zfac = (Zsp[isp]/Zref)^2 + nussp = nuref*zfac + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + C_M_exact[ivpa,ivperp] += Cssp_Maxwellian_inputs(dens,upar,vth,mref, + denssp[isp],uparsp[isp],vthsp[isp],msp[isp], + nussp,vpa,vperp,ivpa,ivperp) + end + end + end + begin_s_r_z_anyv_region() + @views fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!(Fs_M[:,:], + nuref,mref,Zref,msp,Zsp,denssp,uparsp,vthsp, + fkpl_arrays,vperp,vpa,vperp_spectral,vpa_spectral) + if test_numerical_conserving_terms + # 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 + density_conserving_correction!(fkpl_arrays.CC,Fs_M,vpa,vperp,dummy_array) + end + # extract C[Fs,Fs'] result + begin_vperp_vpa_region() + @loop_vperp_vpa ivperp ivpa begin + C_M_num[ivpa,ivperp] = fkpl_arrays.CC[ivpa,ivperp] + end + begin_serial_region() + @serial_region begin + C_M_max, C_M_L2 = print_test_data(C_M_exact,C_M_num,C_M_err,"C_M",vpa,vperp,dummy_array,print_to_screen=print_to_screen) + atol_max = 7.0e-2 + atol_L2 = 6.0e-4 + @test C_M_max < atol_max + @test C_M_L2 < atol_L2 + if !test_numerical_conserving_terms + delta_n = get_density(C_M_num, vpa, vperp) + rtol, atol = 0.0, 1.0e-12 + @test isapprox(delta_n, rtol ; atol=atol) + if print_to_screen + println("delta_n: ", delta_n) + end + elseif test_numerical_conserving_terms + delta_n = get_density(C_M_num, vpa, vperp) + rtol, atol = 0.0, 1.0e-15 + @test isapprox(delta_n, rtol ; atol=atol) + if print_to_screen + println("delta_n: ", delta_n) + end + end + end + end + finalize_comms!() + end + @testset " - test weak-form Rosenbluth potential calculation: direct integration" begin println(" - test weak-form Rosenbluth potential calculation: direct integration") ngrid = 5 # chosen for a quick test -- direct integration is slow! diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index 53b003cac..972b92620 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -2932,12 +2932,31 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(vperp.grid, vpa.grid, pdf[:,:,is,i], xlabel="vperp", ylabel="vpa", c = :deep, interpolation = :cubic) end - outfile = string(run_name_label, "_pdf_vs_vperp_vpa", iz0_string, ir0_string, spec_string[is], ".gif") + outfile = string(run_name_label, "_pdf_vs_vpa_vperp", iz0_string, ir0_string, spec_string[is], ".gif") trygif(anim, outfile, fps=5) @views heatmap(vperp.grid, vpa.grid, pdf[:,:,is,itime_max], xlabel="vperp", ylabel="vpa", c = :deep, interpolation = :cubic) outfile = string(run_name_label, "_pdf_vs_vpa_vperp", ir0_string, iz0_string, spec_string[is], ".pdf") savefig(outfile) + + ivpa = floor(mk_int,vpa.n/2) + 1 + @views plot(vperp.grid, pdf[ivpa,:,is,itime_max], xlabel="vperp", ylabel="f", label="", linewidth=2) + outfile = string(run_name_label, "_pdf_vs_vperp","_ivpa", ivpa, ir0_string, iz0_string, spec_string[is], ".pdf") + savefig(outfile) + ivperp = 1 + @views plot(vpa.grid, pdf[:,ivperp,is,itime_max], xlabel="vpa", ylabel="f", label="", linewidth=2) + outfile = string(run_name_label, "_pdf_vs_vpa", "_ivperp", ivperp, ir0_string, iz0_string, spec_string[is], ".pdf") + savefig(outfile) + + + ivpa = floor(mk_int,vpa.n/2) + 1 + @views plot(vperp.grid, log.(abs.(pdf[ivpa,:,is,itime_max])), xlabel="vperp", ylabel="log|f|", label="", linewidth=2) + outfile = string(run_name_label, "_logpdf_vs_vperp","_ivpa", ivpa, ir0_string, iz0_string, spec_string[is], ".pdf") + savefig(outfile) + ivperp = 1 + @views plot(vpa.grid, log.(abs.(pdf[:,ivperp,is,itime_max])), xlabel="vpa", ylabel="log|f|", label="", linewidth=2) + outfile = string(run_name_label, "_logpdf_vs_vpa", "_ivperp", ivperp, ir0_string, iz0_string, spec_string[is], ".pdf") + savefig(outfile) end elseif vperp.n == 1 for is ∈ 1:n_species