diff --git a/src/input_structs.jl b/src/input_structs.jl index b3cb44238..ed3d20914 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -61,6 +61,7 @@ mutable struct advance_info ionization_collisions::Bool ionization_collisions_1V::Bool ionization_source::Bool + krook_collisions::Bool numerical_dissipation::Bool source_terms::Bool continuity::Bool @@ -305,6 +306,8 @@ mutable struct collisions_input ionization::mk_float # if constant_ionization_rate = true, use an ionization term that is constant in z constant_ionization_rate::Bool + # Coulomb collision rate at the reference density and temperature + krook_collision_frequency_prefactor::mk_float end """ diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl new file mode 100644 index 000000000..a4c40b5ee --- /dev/null +++ b/src/krook_collisions.jl @@ -0,0 +1,129 @@ +""" +""" +module krook_collisions + +using ..constants: epsilon0, proton_charge +using ..looping + +""" +Calculate normalized collision frequency at reference parameters for Coulomb collisions. + +Currently valid only for hydrogenic ions (Z=1) +""" +function setup_krook_collisions(reference_parameters) + Nref = reference_parameters.Nref + Tref = reference_parameters.Tref + mref = reference_parameters.mref + timeref = reference_parameters.timeref + + Nref_per_cm3 = Nref * 1.0e-6 + + # Coulomb logarithm at reference parameters for same-species ion-ion collisions, using + # NRL formulary. Formula given for n in units of cm^-3 and T in units of eV. + logLambda_ii = 23.0 - log(sqrt(2.0*Nref_per_cm3) / Tref^1.5) + + # Collision frequency, using \hat{\nu} from Appendix, p. 277 of Helander "Collisional + # Transport in Magnetized Plasmas" (2002). + T0_Joules = Tref * proton_charge # Tref in Joules + mi_kg = mref # mi in kg + vth_i0 = sqrt(2.0 * T0_Joules / mi_kg) # vth_i at reference parameters in m.s^-1 + nu_ii0_per_s = Nref * proton_charge^4 * logLambda_ii / + (4.0 * π * epsilon0^2 * mi_kg^2 * vth_i0^3) # s^-1 + nu_ii0 = nu_ii0_per_s * timeref + + return nu_ii0 +end + +""" +Add collision operator + +Currently Krook collisions +""" +function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, vperp, vpa, dt) + begin_s_r_z_region() + + if vperp.n > 1 + error("Krook collisions not implemented for 2V case yet") + end + + # Note: do not need 1/sqrt(pi) for the 'Maxwellian' term because the pdf is already + # normalized by sqrt(pi) (see velocity_moments.integrate_over_vspace). + + if moments.evolve_ppar && moments.evolve_upar + # Compared to evolve_upar version, grid is already normalized by vth, and multiply + # through by vth, remembering pdf is already multiplied by vth + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + T = (moments.charged.vth[iz,ir,is])^2 + nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - exp(-vpa.grid[ivpa]^2 - vperp.grid[ivperp]^2)) + end + end + elseif moments.evolve_ppar + # Compared to full-f collision operater, multiply through by vth, remembering pdf + # is already multiplied by vth, and grid is already normalized by vth + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + T = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is])/vth)^2 + - (vperp.grid[ivperp]/vth)^2)) + end + end + elseif moments.evolve_upar + # Compared to evolve_density version, grid is already shifted by upar + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + T = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - 1.0 / vth * exp(-(vpa.grid[ivpa] / vth)^2 + - (vperp.grid[ivperp] / vth)^2)) + end + end + elseif moments.evolve_density + # Compared to full-f collision operater, divide through by density, remembering + # that pdf is already normalized by density + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + T = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - 1.0 / vth + * exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is]) / vth)^2 + - (vperp.grid[ivperp]/vth)^2)) + end + end + else + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + T = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - n / vth + * exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is])/vth)^2 + - (vperp.grid[ivperp]/vth)^2)) + end + end + end + + return nothing +end + +end # krook_collisions diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 73bfd15b8..03c0a557f 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -46,6 +46,7 @@ include("neutral_z_advection.jl") include("neutral_vz_advection.jl") include("charge_exchange.jl") include("ionization.jl") +include("krook_collisions.jl") include("continuity.jl") include("energy_equation.jl") include("force_balance.jl") diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 31f0d61e4..752fb68c5 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -13,6 +13,7 @@ using ..communication using ..coordinates: define_coordinate using ..file_io: io_has_parallel, input_option_error, open_ascii_output_file using ..constants +using ..krook_collisions: setup_krook_collisions using ..finite_differences: fd_check_option using ..input_structs using ..numerical_dissipation: setup_numerical_dissipation @@ -177,6 +178,12 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) collisions.charge_exchange = get(scan_input, "charge_exchange_frequency", 2.0*sqrt(species.charged[1].initial_temperature)) collisions.ionization = get(scan_input, "ionization_frequency", collisions.charge_exchange) collisions.constant_ionization_rate = get(scan_input, "constant_ionization_rate", false) + include_krook_collisions = get(scan_input, "krook_collisions", false) + if include_krook_collisions + collisions.krook_collision_frequency_prefactor = setup_krook_collisions(reference_parameters) + else + collisions.krook_collision_frequency_prefactor = -1.0 + end # parameters related to the time stepping nstep = get(scan_input, "nstep", 5) @@ -948,7 +955,9 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # ionization collision frequency ionization = 0.0 constant_ionization_rate = false - collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate) + krook_collision_frequency_prefactor = -1.0 + collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, + krook_collision_frequency_prefactor) Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength diff --git a/src/time_advance.jl b/src/time_advance.jl index 9d2a12d99..4892dd87e 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -40,6 +40,7 @@ using ..vperp_advection: update_speed_vperp!, vperp_advection! using ..vpa_advection: update_speed_vpa!, vpa_advection! using ..charge_exchange: charge_exchange_collisions_1V!, charge_exchange_collisions_3V! using ..ionization: ionization_collisions_1V!, ionization_collisions_3V!, constant_ionization_source! +using ..krook_collisions: krook_collisions! using ..numerical_dissipation: vpa_boundary_buffer_decay!, vpa_boundary_buffer_diffusion!, vpa_dissipation!, z_dissipation!, r_dissipation!, @@ -404,6 +405,7 @@ function setup_advance_flags(moments, composition, t_input, collisions, num_diss advance_ionization = false advance_ionization_1V = false advance_ionization_source = false + advance_krook_collisions = false advance_numerical_dissipation = false advance_sources = false advance_continuity = false @@ -473,6 +475,9 @@ function setup_advance_flags(moments, composition, t_input, collisions, num_diss if collisions.ionization > 0.0 && collisions.constant_ionization_rate advance_ionization_source = true end + if collisions.krook_collision_frequency_prefactor > 0.0 + advance_krook_collisions = true + end advance_numerical_dissipation = true # if evolving the density, must advance the continuity equation, # in addition to including sources arising from the use of a modified distribution @@ -521,12 +526,12 @@ function setup_advance_flags(moments, composition, t_input, collisions, num_diss advance_neutral_z_advection, advance_neutral_r_advection, advance_neutral_vz_advection, advance_cx, advance_cx_1V, advance_ionization, advance_ionization_1V, - advance_ionization_source, advance_numerical_dissipation, - advance_sources, advance_continuity, advance_force_balance, - advance_energy, advance_neutral_sources, - advance_neutral_continuity, advance_neutral_force_balance, - advance_neutral_energy, rk_coefs, manufactured_solns_test, - r_diffusion, vpa_diffusion, vz_diffusion) + advance_ionization_source, advance_krook_collisions, + advance_numerical_dissipation, advance_sources, + advance_continuity, advance_force_balance, advance_energy, + advance_neutral_sources, advance_neutral_continuity, + advance_neutral_force_balance, advance_neutral_energy, rk_coefs, + manufactured_solns_test, r_diffusion, vpa_diffusion, vz_diffusion) end function setup_dummy_and_buffer_arrays(nr,nz,nvpa,nvperp,nvz,nvr,nvzeta,nspecies_ion,nspecies_neutral) @@ -1039,6 +1044,14 @@ function time_advance_split_operators!(pdf, scratch, t, t_input, vpa, z, advance.ionization_collisions = false end end + if collisions.krook_collision_frequency_prefactor > 0.0 + advance.krook_collisions = true + time_advance_no_splitting!(pdf, scratch, t, t_input, z, vpa, + z_spectral, vpa_spectral, moments, fields, z_advect, vpa_advect, + z_SL, vpa_SL, composition, collisions, sources, num_diss_params, + advance, istep) + advance.krook_collisions = false + end # and add the source terms associated with redefining g = pdf/density or pdf*vth/density # to the kinetic equation if moments.evolve_density || moments.evolve_upar || moments.evolve_ppar @@ -1652,6 +1665,12 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, collisions, dt) end + # Add a for Krook collision operoator for ions + if advance.krook_collisions + krook_collisions!(fvec_out.pdf, fvec_in, moments, composition, collisions, + vperp, vpa, dt) + end + # add numerical dissipation if advance.numerical_dissipation vpa_dissipation!(fvec_out.pdf, fvec_in.pdf, vpa, vpa_spectral, dt, diff --git a/src/utils.jl b/src/utils.jl index fcf549b69..5f9c50636 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -63,6 +63,8 @@ function get_unnormalized_parameters(input::Dict) parameters["CX_rate_coefficient"] = collisions.charge_exchange / Nnorm / timenorm parameters["ionization_rate_coefficient"] = collisions.ionization / Nnorm / timenorm + parameters["coulomb_collision_frequency0"] = + collisions.coulomb_collision_frequency_prefactor / timenorm return parameters end