Krook collision operator for ion-ion Coulomb collisions

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
 
 """

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

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")

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

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,

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