Merge pull request #240 from mabarnes/remove-constant-ionization-source

Remove 'constant ionization source'

moment_kinetics/debug_test/runtests.jl

@@ -8,7 +8,6 @@ function runtests()
         include(joinpath(@__DIR__, "sound_wave_tests.jl"))
         include(joinpath(@__DIR__, "fokker_planck_collisions_tests.jl"))
         include(joinpath(@__DIR__, "wall_bc_tests.jl"))
-        include(joinpath(@__DIR__, "harrisonthompson.jl"))
         #include(joinpath(@__DIR__, "mms_tests.jl"))
         include(joinpath(@__DIR__, "restart_interpolation_tests.jl"))
         include(joinpath(@__DIR__, "recycling_fraction_tests.jl"))

moment_kinetics/src/input_structs.jl

@@ -115,7 +115,6 @@ mutable struct advance_info
     neutral_ionization_collisions::Bool
     ion_ionization_collisions_1V::Bool
     neutral_ionization_collisions_1V::Bool
-    ionization_source::Bool
     krook_collisions_ii::Bool
     mxwl_diff_collisions_ii::Bool
     mxwl_diff_collisions_nn::Bool
@@ -458,8 +457,6 @@ struct collisions_input
     ionization_electron::mk_float
     # ionization energy cost
     ionization_energy::mk_float
-    # if constant_ionization_rate = true, use an ionization term that is constant in z
-    constant_ionization_rate::Bool
     # electron-ion collision frequency
     nu_ei::mk_float
     # struct of parameters for the Krook operator

moment_kinetics/src/ionization.jl

@@ -6,66 +6,10 @@ export ion_ionization_collisions_1V!
 export neutral_ionization_collisions_1V!
 export ion_ionization_collisions_3V!
 export neutral_ionization_collisions_3V!
-export constant_ionization_source!
 
 using ..interpolation: interpolate_to_grid_vpa!
 using ..looping
 
-"""
-"""
-
-function constant_ionization_source!(f_out, fvec_in, vpa, vperp, z, r, moments,
-                                     composition, collisions, dt)
-    @boundscheck vpa.n == size(f_out,1) || throw(BoundsError(f_out))
-    @boundscheck vperp.n == size(f_out,2) || throw(BoundsError(f_out))
-    @boundscheck z.n == size(f_out,3) || throw(BoundsError(f_out))
-    @boundscheck r.n == size(f_out,4) || throw(BoundsError(f_out))
-    @boundscheck composition.n_ion_species == size(f_out,5) || throw(BoundsError(f_out))
-
-    begin_s_r_z_region()
-
-    # Oddly the test in test/harrisonthompson.jl matches the analitical
-    # solution (which assumes width=0.0) better with width=0.5 than with,
-    # e.g., width=0.15. Possibly narrower widths would require more vpa
-    # resolution, which then causes crashes due to overshoots giving
-    # negative f??
-    width = 0.5
-    vperpwidth = 0.5
-    rwidth = 0.5
-    if vperp.n > 1
-        vperpprefac = 1.0/vperpwidth^2
-    else
-        vperpprefac = 1.0
-    end
-    # loop below relies on vperp[1] = 0 when vperp.n = 1
-    @loop_s_r is ir begin
-        rfac = exp( - (r.grid[ir]/rwidth)^2)
-
-        @loop_z iz begin
-            if moments.evolve_ppar && moments.evolve_upar
-                @. vpa.scratch = vpa.grid / moments.vth[iz] + fvec_in.upar[iz]
-                prefactor = moments.vth[iz] / fvec_in.dens[iz]
-            elseif moments.evolve_ppar
-                @. vpa.scratch = vpa.grid / moments.vth[iz]
-                prefactor = moments.vth[iz] / fvec_in.dens[iz]
-            elseif moments.evolve_upar
-                @. vpa.scratch = vpa.grid + fvec_in.upar[iz]
-                prefactor = 1.0 / fvec_in.dens[iz]
-            elseif moments.evolve_density
-                @. vpa.scratch = vpa.grid
-                prefactor = 1.0 / fvec_in.dens[iz]
-            else
-                @. vpa.scratch = vpa.grid
-                prefactor = 1.0
-            end
-            @loop_vperp_vpa ivperp ivpa begin
-                vperpfac = vperpprefac*exp( - (vperp.grid[ivperp]/vperpwidth)^2) 
-                f_out[ivpa,ivperp,iz,ir,is] += dt*rfac*vperpfac*collisions.ionization/width*prefactor*exp(-(vpa.scratch[ivpa]/width)^2)
-            end
-        end
-    end
-end 
-
 function ion_ionization_collisions_1V!(f_out, fvec_in, vz, vpa, vperp, z, r, vz_spectral,
                                        moments, composition, collisions, dt)
     # This routine assumes a 1D model with:

moment_kinetics/src/moment_kinetics_input.jl

@@ -190,7 +190,6 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     ionization = get(scan_input, "ionization_frequency", charge_exchange)
     ionization_electron = get(scan_input, "electron_ionization_frequency", ionization)
     ionization_energy = get(scan_input, "ionization_energy", 0.0)
-    constant_ionization_rate = get(scan_input, "constant_ionization_rate", false)
     nu_ei = get(scan_input, "nu_ei", 0.0)
     # set up krook collision inputs
     krook_input = setup_krook_collisions_input(scan_input, reference_params)
@@ -202,10 +201,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # for the collisions outputs itself a struct of the type of collision, which
     # is a substruct of the overall collisions_input struct.
     collisions = collisions_input(charge_exchange, charge_exchange_electron, ionization,
-                                  ionization_electron, ionization_energy,
-                                  constant_ionization_rate, 
-                                  nu_ei, krook_input,
-                                  fkpl_input, mxwl_diff_input)
+                                  ionization_electron, ionization_energy, nu_ei,
+                                  krook_input, fkpl_input, mxwl_diff_input)
 
     # parameters related to the time stepping
     timestepping_section = set_defaults_and_check_section!(

moment_kinetics/src/time_advance.jl

@@ -48,8 +48,7 @@ using ..charge_exchange: ion_charge_exchange_collisions_1V!,
                          neutral_charge_exchange_collisions_3V!
 using ..electron_kinetic_equation: update_electron_pdf!, implicit_electron_advance!
 using ..ionization: ion_ionization_collisions_1V!, neutral_ionization_collisions_1V!,
-                    ion_ionization_collisions_3V!, neutral_ionization_collisions_3V!,
-                    constant_ionization_source!
+                    ion_ionization_collisions_3V!, neutral_ionization_collisions_3V!
 using ..krook_collisions: krook_collisions!
 using ..maxwell_diffusion: ion_vpa_maxwell_diffusion!, neutral_vz_maxwell_diffusion!
 using ..external_sources
@@ -1106,7 +1105,6 @@ function setup_advance_flags(moments, composition, t_params, collisions,
     advance_neutral_ionization = false
     advance_ion_ionization_1V = false
     advance_neutral_ionization_1V = false
-    advance_ionization_source = false
     advance_krook_collisions_ii = false
     advance_maxwell_diffusion_ii = false
     advance_maxwell_diffusion_nn = false
@@ -1188,10 +1186,6 @@ function setup_advance_flags(moments, composition, t_params, collisions,
                 end
             end
         end
-        # exception for the case where ions are evolved alone but sourced by ionization
-        if collisions.ionization > 0.0 && collisions.constant_ionization_rate && !t_params.implicit_ion_advance
-            advance_ionization_source = true
-        end
         # set flags for krook and maxwell diffusion collisions, and negative coefficient
         # in both cases (as usual) will mean not employing that operator (flag remains false)
         if collisions.krook.nuii0 > 0.0
@@ -1283,8 +1277,7 @@ function setup_advance_flags(moments, composition, t_params, collisions,
                         advance_neutral_vz_advection, advance_ion_cx, advance_neutral_cx,
                         advance_ion_cx_1V, advance_neutral_cx_1V, advance_ion_ionization,
                         advance_neutral_ionization, advance_ion_ionization_1V,
-                        advance_neutral_ionization_1V, advance_ionization_source,
-                        advance_krook_collisions_ii,
+                        advance_neutral_ionization_1V, advance_krook_collisions_ii,
                         advance_maxwell_diffusion_ii, advance_maxwell_diffusion_nn,
                         explicit_weakform_fp_collisions,
                         advance_external_source, advance_ion_numerical_dissipation,
@@ -1320,7 +1313,6 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions
     advance_neutral_ionization = false
     advance_ion_ionization_1V = false
     advance_neutral_ionization_1V = false
-    advance_ionization_source = false
     advance_krook_collisions_ii = false
     advance_maxwell_diffusion_ii = false
     advance_maxwell_diffusion_nn = false
@@ -1376,7 +1368,6 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions
                       * "vpa.n=$(vpa.n), vz.n=$(vz.n)")
             end
         end
-        advance_ionization_source = collisions.ionization > 0.0 && collisions.constant_ionization_rate
         advance_krook_collisions_ii = collisions.krook.nuii0 > 0.0
         advance_external_source = external_source_settings.ion.active
         advance_ion_numerical_dissipation = true
@@ -1418,8 +1409,7 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions
                         advance_neutral_vz_advection, advance_ion_cx, advance_neutral_cx,
                         advance_ion_cx_1V, advance_neutral_cx_1V, advance_ion_ionization,
                         advance_neutral_ionization, advance_ion_ionization_1V,
-                        advance_neutral_ionization_1V,
-                        advance_ionization_source, advance_krook_collisions_ii,
+                        advance_neutral_ionization_1V, advance_krook_collisions_ii,
                         advance_maxwell_diffusion_ii, advance_maxwell_diffusion_nn,
                         explicit_weakform_fp_collisions,
                         advance_external_source, advance_ion_numerical_dissipation,
@@ -3342,10 +3332,6 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments,
         neutral_ionization_collisions_3V!(fvec_out.pdf_neutral, fvec_in, composition, vz,
                                           vr, vzeta, vpa, vperp, z, r, collisions, dt)
     end
-    if advance.ionization_source
-        constant_ionization_source!(fvec_out.pdf, fvec_in, vpa, vperp, z, r, moments,
-                                    composition, collisions, dt)
-    end
 
     # Add Krook collision operator for ions
     if advance.krook_collisions_ii