From d17e8bf01f632b9d416d6a16fa6c5382fbb8a379 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 8 Aug 2024 14:12:01 +0100 Subject: [PATCH] Remove 'constant ionization source' This functionality has been replaced with the more flexible 'external sources' feature, so can now be removed. Also removes the 'Harrison-Thompson' debug check as this now includes no features that are not in the 'recycling fraction' check. --- .../debug_test/harrisonthompson.jl | 23 ------ .../debug_test/harrisonthompson_inputs.jl | 76 ------------------- moment_kinetics/debug_test/runtests.jl | 1 - moment_kinetics/src/input_structs.jl | 3 - moment_kinetics/src/ionization.jl | 56 -------------- moment_kinetics/src/moment_kinetics_input.jl | 7 +- moment_kinetics/src/time_advance.jl | 20 +---- 7 files changed, 5 insertions(+), 181 deletions(-) delete mode 100644 moment_kinetics/debug_test/harrisonthompson.jl delete mode 100644 moment_kinetics/debug_test/harrisonthompson_inputs.jl diff --git a/moment_kinetics/debug_test/harrisonthompson.jl b/moment_kinetics/debug_test/harrisonthompson.jl deleted file mode 100644 index b310a7d11..000000000 --- a/moment_kinetics/debug_test/harrisonthompson.jl +++ /dev/null @@ -1,23 +0,0 @@ -module HarrisonThompsonDebug -# Test case with constant-in-space, delta-function-in-vpa source against -# analytic solution from [E. R. Harrison and W. B. Thompson. The low pressure -# plane symmetric discharge. Proc. Phys. Soc., 74:145, 1959] - -include("setup.jl") - -# Create a temporary directory for test output -test_output_directory = get_MPI_tempdir() -mkpath(test_output_directory) - -# Input parameters for the test -include("harrisonthompson_inputs.jl") - -# Defines the test functions, using variables defined in the *_inputs.jl file -include("runtest_template.jl") - -end # HarrisonThompsonDebug - - -using .HarrisonThompsonDebug - -HarrisonThompsonDebug.runtests() diff --git a/moment_kinetics/debug_test/harrisonthompson_inputs.jl b/moment_kinetics/debug_test/harrisonthompson_inputs.jl deleted file mode 100644 index 6b12546a8..000000000 --- a/moment_kinetics/debug_test/harrisonthompson_inputs.jl +++ /dev/null @@ -1,76 +0,0 @@ -test_type = "Harrison-Thompson wall boundary condition" - -# default inputs for tests -test_input_finite_difference = Dict("n_ion_species" => 2, - "n_neutral_species" => 2, - "boltzmann_electron_response" => true, - "run_name" => "finite_difference", - "base_directory" => test_output_directory, - "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" => "gaussian", - "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, - "vpa_IC_option1" => "gaussian", - "vpa_IC_density_amplitude1" => 1.0, - "vpa_IC_density_phase1" => 0.0, - "vpa_IC_upar_amplitude1" => 0.0, - "vpa_IC_upar_phase1" => 0.0, - "vpa_IC_temperature_amplitude1" => 0.0, - "vpa_IC_temperature_phase1" => 0.0, - "charge_exchange_frequency" => 0.0, - #"ionization_frequency" => 0.25, - "ionization_frequency" => 0.688, - "constant_ionization_rate" => true, - "timestepping" => Dict{String,Any}("nstep" => 3, - "dt" => 0.0005, - "nwrite" => 2, - "type" => "SSPRK2", - "split_operators" => false), - "z_ngrid" => 4, - "z_nelement" => 1, - "z_bc" => "periodic", - "z_discretization" => "finite_difference", - "vpa_ngrid" => 4, - "vpa_nelement" => 1, - "vpa_L" => 8.0, - "vpa_bc" => "periodic", - "vpa_discretization" => "finite_difference", - "vperp_ngrid" => 1, - "vperp_nelement" => 1, - "vz_ngrid" => 4, - "vz_nelement" => 1, - "vz_L" => 8.0, - "vz_bc" => "periodic", - "vz_discretization" => "finite_difference", - "vzeta_ngrid" => 1, - "vzeta_nelement" => 1, - "vr_ngrid" => 1, - "vr_nelement" => 1) - -test_input_chebyshev = merge(test_input_finite_difference, - Dict("run_name" => "chebyshev_pseudospectral", - "z_discretization" => "chebyshev_pseudospectral", - "z_ngrid" => 3, - "z_nelement" => 2, - "vpa_discretization" => "chebyshev_pseudospectral", - "vpa_ngrid" => 3, - "vpa_nelement" => 2, - "vz_discretization" => "chebyshev_pseudospectral", - "vz_ngrid" => 3, - "vz_nelement" => 2)) - -test_input_list = [ - #test_input_finite_difference, - test_input_chebyshev, - ] diff --git a/moment_kinetics/debug_test/runtests.jl b/moment_kinetics/debug_test/runtests.jl index 86d3056ac..6e16a0236 100644 --- a/moment_kinetics/debug_test/runtests.jl +++ b/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")) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 57f692c62..30bc39cb2 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/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 diff --git a/moment_kinetics/src/ionization.jl b/moment_kinetics/src/ionization.jl index 9ea348da0..e1becf468 100644 --- a/moment_kinetics/src/ionization.jl +++ b/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: diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 57c1c8205..12e0a879f 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/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!( diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index db4464985..774028179 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/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