From 28905e065bbf24f5bba17b6b4ae15c4fbe1fe4f6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 28 Mar 2024 12:19:39 +0000 Subject: [PATCH 01/26] Port changes to MPI.Init() call. --- moment_kinetics/src/communication.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/communication.jl b/moment_kinetics/src/communication.jl index 072d7b594..a8e4e06d3 100644 --- a/moment_kinetics/src/communication.jl +++ b/moment_kinetics/src/communication.jl @@ -117,8 +117,10 @@ const global_Win_store = Vector{MPI.Win}(undef, 0) """ """ function __init__() - MPI.Init() - + if !MPI.Initialized() + MPI.Init() + end + comm_world.val = MPI.COMM_WORLD.val global_rank[] = MPI.Comm_rank(comm_world) From 21e8c2158ac6a9a2bc82bd78e7f67469323bd169 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 28 Mar 2024 12:20:07 +0000 Subject: [PATCH 02/26] Add element boundaries to coord struct for use in Lagrange interpolation. --- moment_kinetics/src/coordinates.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl index 936695a10..557f5f47b 100644 --- a/moment_kinetics/src/coordinates.jl +++ b/moment_kinetics/src/coordinates.jl @@ -98,6 +98,8 @@ struct coordinate element_shift::Array{mk_float,1} # option used to set up element spacing element_spacing_option::String + # list of element boundaries + element_boundaries::Array{mk_float,1} end """ @@ -166,7 +168,8 @@ function define_coordinate(input, parallel_io::Bool=false; init_YY::Bool=true) 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_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) + 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() From 2478e3edd48c888790cd3274f9f7364c085b9daf Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 28 Mar 2024 12:21:54 +0000 Subject: [PATCH 03/26] Add facility to interpolate and extrapolate pdf data from one grid to another, assuming that the physical meaning of the coordinate values change, but not the normalised values themselves. We assume that the pdf goes to zero by the end of the "natural" normalised grid, and thus can extrapolate the pdf with zeros into the region beyond which data exists. --- moment_kinetics/src/fokker_planck_calculus.jl | 94 +++++++++++++++++++ moment_kinetics/test/fokker_planck_tests.jl | 75 +++++++++++++++ 2 files changed, 169 insertions(+) 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/test/fokker_planck_tests.jl b/moment_kinetics/test/fokker_planck_tests.jl index 54266f518..bf2209056 100644 --- a/moment_kinetics/test/fokker_planck_tests.jl +++ b/moment_kinetics/test/fokker_planck_tests.jl @@ -22,6 +22,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 +71,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") From 77cd7511ea290d859d09266eb18c63916592943b Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 10 Apr 2024 12:15:57 +0100 Subject: [PATCH 04/26] Correct normalised collision frequency documentation to normalise correctly to reference timescale. --- moment_kinetics/src/fokker_planck.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 75e38867d..8a4763ac0 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -227,7 +227,7 @@ The result is stored in the array `fkpl_arrays.CC`. The normalised collision frequency is defined by ```math -\\nu_{ss'} = \\frac{\\gamma_{ss'} n_\\mathrm{ref}}{2 m_s^2 c_\\mathrm{ref}^3} +\\tilde{\\nu}_{ss'} = \\frac{L_\\mathrm{ref}}{c_\\mathrm{ref}}\\frac{\\gamma_{ss'} n_\\mathrm{ref}}{2 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\$. From 6be618707dd99c7413018ed29db275a3c56aa6a9 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 10 Apr 2024 14:31:46 +0100 Subject: [PATCH 05/26] Add special Fokker-Planck collision operator function to mock up fast or beam ion collisions with the existing self collision operator and a Maxwellian (slow) ion and electron species. The special function is used because i) electrons are not available in the branch, ii) we do not wish the slow ions or the electrons to evolve. Noting this, we can construct the function so that only two matrix assembly steps are carried out, and only a single set of elliptic solves are needed (for the self collision operator). This means that this test operator is only a factor of two slower than the existing self-collision operator. Input parameters for the unevolved species are hard coded, but should be later made input parameters in the fokker planck input TOML namelist. Changes towards this are made, but currently they are incomplete. --- ...fokker-planck-relaxation-slowing-down.toml | 63 ++++++ moment_kinetics/src/fokker_planck.jl | 183 ++++++++++++++++++ moment_kinetics/src/input_structs.jl | 11 ++ moment_kinetics/src/moment_kinetics_input.jl | 3 +- moment_kinetics/src/time_advance.jl | 11 +- 5 files changed, 268 insertions(+), 3 deletions(-) create mode 100644 examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml 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..918d943de --- /dev/null +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml @@ -0,0 +1,63 @@ +# 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 +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.0 +slowing_down_test = true +nstep = 1000 +dt = 1.0e-5 +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 +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 + diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 8a4763ac0..3df433690 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 @@ -50,6 +51,7 @@ using ..communication using ..velocity_moments: integrate_over_vspace using ..velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_qpar, get_pressure, get_rmom using ..looping +using ..input_structs: fkpl_collisions_input using ..fokker_planck_calculus: init_Rosenbluth_potential_integration_weights! using ..fokker_planck_calculus: init_Rosenbluth_potential_boundary_integration_weights! using ..fokker_planck_calculus: allocate_boundary_integration_weights @@ -73,6 +75,12 @@ using ..fokker_planck_test: F_Maxwellian, dFdvpa_Maxwellian, dFdvperp_Maxwellian # where the potentials are computed by elliptic solve ######################################################## +function setup_fkpl_collisions_input(toml_input::Dict) + input_section = get(toml_input, "fokker_planck_collisions", Dict{String,Any}()) + input = Dict(Symbol(k)=>v for (k,v) in input_section) + return fkpl_collisions_input(; input...) +end + """ function that initialises the arrays needed for Fokker Planck collisions using numerical integration to compute the Rosenbluth potentials only @@ -157,6 +165,95 @@ 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)) + + # masses and collision frequencies + mref = 1.0 + nuref = collisions.nuii # generalise! + #msp = Array{mk_float,1}(undef,2) + #densp = Array{mk_float,1}(undef,2) + #uparsp = Array{mk_float,1}(undef,2) + #vthsp = Array{mk_float,1}(undef,2) + #msp[1], msp[2] = 0.25, 0.25/1836.0 + #uparsp[1], uparsp[2] = 0.0, 0.0 + msp = [0.25, 0.25/1836.0] + densp = [1.0, 1.0] + uparsp = [0.0, 0.0] + vthsp = [sqrt(0.01/msp[1]), sqrt(0.01/msp[2])] + + # N.B. parallelisation using special 'anyv' region + begin_s_r_z_anyv_region() + @loop_s_r_z is ir iz begin + # first argument is Fs, and second argument is Fs' in C[Fs,Fs'] + @views fokker_planck_collision_operator_weak_form!( + pdf_in[:,:,iz,ir,is], pdf_in[:,:,iz,ir,is], mref, mref, nuref, 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 ad-hoc conserving corrections + conserving_corrections!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, + fkpl_arrays.S_dummy) + + # advance this part of s,r,z with the resulting 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 + if diagnose_entropy_production + # assign dummy array + lnfC = fkpl_arrays.rhsvpavperp + @loop_vperp_vpa ivperp ivpa begin + lnfC[ivpa,ivperp] = log(abs(pdf_in[ivpa,ivperp,iz,ir,is]) + 1.0e-15)*CC[ivpa,ivperp] + end + begin_anyv_region() + @anyv_serial_region begin + dSdt[iz,ir,is] = -get_density(lnfC,vpa,vperp) + end + end + + # 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,msp,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) + # advance this part of s,r,z with the resulting 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 """ @@ -341,6 +438,92 @@ 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, + msp::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 + for isp in 1:nsp + dens = densp[isp] + upar = uparsp[isp] + vth = vthsp[isp] + @loop_vperp_vpa ivperp ivpa begin + d2Gdvpa2[ivpa,ivperp] += d2Gdvpa2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + d2Gdvperp2[ivpa,ivperp] += d2Gdvperp2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + d2Gdvperpdvpa[ivpa,ivperp] += d2Gdvperpdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + dHdvpa[ivpa,ivperp] += (ms/msp[isp])*dHdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + dHdvperp[ivpa,ivperp] += (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 diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 8e8e67bbe..462524027 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -323,6 +323,17 @@ mutable struct collisions_input # nu_{ss'} = gamma_{ss'} n_{ref} / 2 (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 nuii::mk_float + # option to determine which FP collision operator is used + slowing_down_test::Bool +end + +""" +""" +Base.@kwdef struct fkpl_collisions_input + # ion-ion self collision frequency + # nu_{ss'} = gamma_{ss'} n_{ref} / 2 (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 + nuii::mk_float end """ diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 691571fca..dc4e96fa6 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -192,6 +192,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) end # set the Fokker-Planck collision frequency collisions.nuii = get(scan_input, "nuii", 0.0) + collisions.slowing_down_test = get(scan_input, "slowing_down_test", false) # parameters related to the time stepping nstep = get(scan_input, "nstep", 5) @@ -1014,7 +1015,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) krook_collision_frequency_prefactor = -1.0 nuii = 0.0 collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, - krook_collision_frequency_prefactor,"none", nuii) + krook_collision_frequency_prefactor,"none", nuii, false) return z, r, vpa, vperp, gyrophase, vz, vr, vzeta, species, composition, drive, evolve_moments, collisions end diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index acfcb6a43..407f5ac3c 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -55,6 +55,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 ..manufactured_solns: manufactured_sources using ..advection: advection_info using ..utils: to_minutes @@ -1764,9 +1765,15 @@ 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.charged.dSdt,composition,collisions,dt, - fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, + if collisions.slowing_down_test + explicit_fp_collisions_weak_form_Maxwellian_cross_species!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt,composition,collisions,dt, + fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral; diagnose_entropy_production = update_entropy_diagnostic) + else + explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt,composition,collisions,dt, + fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, + diagnose_entropy_production = update_entropy_diagnostic) + end end # End of advance for distribution function From 658fd1c33b5c484ff2ed10a8c53cda07f2caa4a6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 10 Apr 2024 15:24:08 +0100 Subject: [PATCH 06/26] Add numerical density conserving term for the cross collision operator. --- moment_kinetics/src/fokker_planck.jl | 35 +++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 3df433690..60efb421d 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -242,7 +242,10 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ 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) - # advance this part of s,r,z with the resulting C[Fs,Fs] + # 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 @@ -619,6 +622,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 From 62503805d7b27bb8e249ba4a48c27f87f2eb5073 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 10 Apr 2024 15:41:02 +0100 Subject: [PATCH 07/26] New initial condition option to mock up an isotropic beam in vpa vperp space. --- ...r-planck-relaxation-slowing-down-beam.toml | 64 +++++++++++++++++++ moment_kinetics/src/initial_conditions.jl | 12 ++++ 2 files changed, 76 insertions(+) create mode 100644 examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml 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..ac35c8a81 --- /dev/null +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml @@ -0,0 +1,64 @@ +# 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" +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.0 +slowing_down_test = true +nstep = 1000 +dt = 1.0e-5 +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 +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 + diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index d0500c986..48c7af3b0 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -690,6 +690,18 @@ function init_charged_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 = 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2) # birth speed of beam + vslow = 0.2*v0 # spread of the beam in speed + v4norm = (vslow^2) * (v0^2) + @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 end return nothing end From 3ced39e99aadc8974d3b6910195be21620a4590a Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 Apr 2024 15:18:22 +0100 Subject: [PATCH 08/26] Bugfix for extracting input collision frequency following merge. --- moment_kinetics/src/fokker_planck.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index f5adefd80..68a6d4001 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -232,7 +232,7 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ # masses and collision frequencies mref = 1.0 - nuref = collisions.nuii # generalise! + nuref = collisions.fkpl.nuii # generalise! #msp = Array{mk_float,1}(undef,2) #densp = Array{mk_float,1}(undef,2) #uparsp = Array{mk_float,1}(undef,2) From 093e93f04c52d6d51579cc8ad62d4cdc12d3355c Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 Apr 2024 15:18:54 +0100 Subject: [PATCH 09/26] Update example input files following merge. --- .../fokker-planck-relaxation-slowing-down-beam.toml | 9 ++++++--- .../fokker-planck-relaxation-slowing-down.toml | 9 ++++++--- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml index ac35c8a81..7128f4173 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml @@ -30,9 +30,6 @@ vpa_IC_option1 = "isotropic-beam" charge_exchange_frequency = 0.0 ionization_frequency = 0.0 constant_ionization_rate = false -# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency -nuii = 1.0 -slowing_down_test = true nstep = 1000 dt = 1.0e-5 nwrite = 50 @@ -62,3 +59,9 @@ 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 diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml index 918d943de..58b8d66e1 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down.toml @@ -29,9 +29,6 @@ z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 constant_ionization_rate = false -# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency -nuii = 1.0 -slowing_down_test = true nstep = 1000 dt = 1.0e-5 nwrite = 50 @@ -61,3 +58,9 @@ 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 From d999c7437fa8e52e6e1d8eb8f576b41dda0bc317 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 Apr 2024 16:57:11 +0100 Subject: [PATCH 10/26] Additional velocity space initial condition for looking at anisotropic beam initial conditions. --- moment_kinetics/src/initial_conditions.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 48c7af3b0..a0849f17f 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -702,6 +702,19 @@ function init_charged_pdf_over_density!(pdf, spec, composition, vpa, vperp, z, 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 = 0.25*0.5*abs(vpa.L) # centre of beam in vpa + vperp0 = 0.5*abs(vperp.L) # centre of beam in vperp + 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 From 24abc2248a18f5ca7a52d6d7f911b6f81d6d88f7 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 Apr 2024 17:19:53 +0100 Subject: [PATCH 11/26] Include Zsprime^2 factor from gamma_ssprime in the cross collision operator, using the linearity of the collision operator in Fsp to sum the Rosenbluth potential contributions into a single array that can be passed to the assembly routine once. The collsion operator should now have the correct O(1) collision frequencies for the cross collisions in normalised units. --- moment_kinetics/src/fokker_planck.jl | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 68a6d4001..b73b4fb7f 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -230,8 +230,9 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ @boundscheck r.n == size(dSdt,2) || throw(BoundsError(dSdt)) @boundscheck n_ion_species == size(dSdt,3) || throw(BoundsError(dSdt)) - # masses and collision frequencies + # masses charge numbers and collision frequencies mref = 1.0 + Zref = 1.0 nuref = collisions.fkpl.nuii # generalise! #msp = Array{mk_float,1}(undef,2) #densp = Array{mk_float,1}(undef,2) @@ -240,6 +241,7 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ #msp[1], msp[2] = 0.25, 0.25/1836.0 #uparsp[1], uparsp[2] = 0.0, 0.0 msp = [0.25, 0.25/1836.0] + Zsp = [0.5, 0.5] densp = [1.0, 1.0] uparsp = [0.0, 0.0] vthsp = [sqrt(0.01/msp[1]), sqrt(0.01/msp[2])] @@ -277,7 +279,7 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ # 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,msp,densp,uparsp,vthsp, + nuref,mref,Zref,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) @@ -492,8 +494,9 @@ 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, - msp::Array{mk_float,1},densp::Array{mk_float,1}, + 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) @@ -530,17 +533,22 @@ function fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!(ffs_in, dHdvpa[ivpa,ivperp] = 0.0 dHdvperp[ivpa,ivperp] = 0.0 end - # sum the contributions from the potentials + # 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]/Zz)^2 # factor from gamma_ss' @loop_vperp_vpa ivperp ivpa begin - d2Gdvpa2[ivpa,ivperp] += d2Gdvpa2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) - d2Gdvperp2[ivpa,ivperp] += d2Gdvperp2_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) - d2Gdvperpdvpa[ivpa,ivperp] += d2Gdvperpdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) - dHdvpa[ivpa,ivperp] += (ms/msp[isp])*dHdvpa_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) - dHdvperp[ivpa,ivperp] += (ms/msp[isp])*dHdvperp_Maxwellian(dens,upar,vth,vpa,vperp,ivpa,ivperp) + 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 From 202b6fae1491b7e76fa7b0093ee811e1d19fc7c5 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 Apr 2024 17:25:57 +0100 Subject: [PATCH 12/26] Remove excess :: --- moment_kinetics/src/fokker_planck.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index b73b4fb7f..1b03fca3f 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -495,7 +495,7 @@ 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}, + 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, From 4ef9014d5a790a312cb9567f95d22bfcc82225fd Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 Apr 2024 17:34:08 +0100 Subject: [PATCH 13/26] Bugfix for ZZ. --- moment_kinetics/src/fokker_planck.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 1b03fca3f..711c65fad 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -542,7 +542,7 @@ function fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!(ffs_in, dens = densp[isp] upar = uparsp[isp] vth = vthsp[isp] - ZZ = (Zsp[isp]/Zz)^2 # factor from gamma_ss' + 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) From 5961076653b9e2522ab2a93db8a193ce4a79d096 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 15 Apr 2024 10:12:54 +0100 Subject: [PATCH 14/26] Extract self-collisions code from cross-collision operator to improved testability of formulation. We envisage that there should be two operators - the self collision operator already implemented and a cross collision operator that sums over all the Rosenbluth potentials of other species in the system. This is because self collisions should include the numerical conserving terms for n, u and T whereas cross collisions can only be forced to conserve n. --- moment_kinetics/src/fokker_planck.jl | 28 ---------------------------- moment_kinetics/src/time_advance.jl | 13 +++++++------ 2 files changed, 7 insertions(+), 34 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 711c65fad..c8405ffb7 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -249,34 +249,6 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ # N.B. parallelisation using special 'anyv' region begin_s_r_z_anyv_region() @loop_s_r_z is ir iz begin - # first argument is Fs, and second argument is Fs' in C[Fs,Fs'] - @views fokker_planck_collision_operator_weak_form!( - pdf_in[:,:,iz,ir,is], pdf_in[:,:,iz,ir,is], mref, mref, nuref, 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 ad-hoc conserving corrections - conserving_corrections!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, - fkpl_arrays.S_dummy) - - # advance this part of s,r,z with the resulting 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 - if diagnose_entropy_production - # assign dummy array - lnfC = fkpl_arrays.rhsvpavperp - @loop_vperp_vpa ivperp ivpa begin - lnfC[ivpa,ivperp] = log(abs(pdf_in[ivpa,ivperp,iz,ir,is]) + 1.0e-15)*CC[ivpa,ivperp] - end - begin_anyv_region() - @anyv_serial_region begin - dSdt[iz,ir,is] = -get_density(lnfC,vpa,vperp) - end - end - # 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,Zref,msp,Zsp,densp,uparsp,vthsp, diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index fee569472..d6c5ada55 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -1765,14 +1765,15 @@ 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) + # self collisions for each species + explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt,composition, + collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, + diagnose_entropy_production = update_entropy_diagnostic) if collisions.fkpl.slowing_down_test - explicit_fp_collisions_weak_form_Maxwellian_cross_species!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt,composition,collisions,dt, - fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral; + # include cross-collsions with fixed Maxwellian backgrounds + explicit_fp_collisions_weak_form_Maxwellian_cross_species!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt, + composition,collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral, diagnose_entropy_production = update_entropy_diagnostic) - else - explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt,composition,collisions,dt, - fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, - diagnose_entropy_production = update_entropy_diagnostic) end end From d38b90bda9b2d83c22ef87f1c3028b585a38c1a2 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 15 Apr 2024 14:04:24 +0100 Subject: [PATCH 15/26] Add test which explicitly tests the new function for slowing-down collisions on fixed Maxwellian background distribution functions. --- moment_kinetics/test/fokker_planck_tests.jl | 99 ++++++++++++++++++++- 1 file changed, 97 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/test/fokker_planck_tests.jl b/moment_kinetics/test/fokker_planck_tests.jl index bf2209056..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 @@ -217,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) @@ -365,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) @@ -515,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! From a63f86003b9ea5bdc303503e7d7f6f5c3cdee754 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 15 Apr 2024 15:29:05 +0100 Subject: [PATCH 16/26] Add additional options to the Fokker Planck operator for physics tests with slowing down problems. --- moment_kinetics/src/fokker_planck.jl | 31 ++++++++++++++++------------ moment_kinetics/src/input_structs.jl | 16 +++++++++++++- moment_kinetics/src/time_advance.jl | 10 +++++---- 3 files changed, 39 insertions(+), 18 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index c8405ffb7..7bba14338 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -84,8 +84,16 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) # get reference collision frequency (note factor of 1/4 due to definition choices) nuii_fkpl_default = 0.25*get_reference_collision_frequency(reference_params) # read the input toml and specify a sensible default - default_dict = Dict{String,Any}("use_fokker_planck" => false, "nuii" => -1.0, - "frequency_option" => "manual", "slowing_down_test" => false ) + default_dict = Dict{String,Any}("use_fokker_planck" => false, + "nuii" => -1.0, + "frequency_option" => "manual", + "self_collisions"=> true, + "slowing_down_test" => false, + "sd_density" => 1.0, + "sd_temp" => 0.01, + "sd_q" => 0.5, + "sd_mi" => 0.25, + "sd_me" => 0.25/1836.0) input_section = get(toml_input, "fokker_planck_collisions", default_dict) # ensure defaults are carried over from the default dict if only a partial dict is given as an input # as the default Dict here is entirely ignored if the namelist is present @@ -230,21 +238,18 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ @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 Zref = 1.0 - nuref = collisions.fkpl.nuii # generalise! - #msp = Array{mk_float,1}(undef,2) - #densp = Array{mk_float,1}(undef,2) - #uparsp = Array{mk_float,1}(undef,2) - #vthsp = Array{mk_float,1}(undef,2) - #msp[1], msp[2] = 0.25, 0.25/1836.0 - #uparsp[1], uparsp[2] = 0.0, 0.0 - msp = [0.25, 0.25/1836.0] - Zsp = [0.5, 0.5] - densp = [1.0, 1.0] + nuref = fkin.nuii # generalise! + msp = [fkin.sd_mi, fkin.sd_me] + Zsp = [fkin.sd_q, fkin.sd_q] + # assume here that ne = sum_i n_i and that initial condition + # for beam ions has unit density + densp = [fkin.sd_density, fkin.sd_density+1.0] uparsp = [0.0, 0.0] - vthsp = [sqrt(0.01/msp[1]), sqrt(0.01/msp[2])] + 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() diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index a0c473c7e..0b25dcbce 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -317,15 +317,29 @@ 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'} = (L/c_{ref}) * gamma_{ss'} n_{ref} / 2 (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 nuii::mk_float - # option to determine which FP collision operator is used + # 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 + sd_q::mk_float + # ion mass with respect to reference + sd_mi::mk_float + # electron mass with respect to reference + sd_me::mk_float end """ diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index d6c5ada55..7864211fe 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -1765,10 +1765,12 @@ 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) - # self collisions for each species - explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.charged.dSdt,composition, - collisions,dt,fp_arrays,r,z,vperp,vpa,vperp_spectral,vpa_spectral,scratch_dummy, - diagnose_entropy_production = update_entropy_diagnostic) + if collisions.fkpl.self_collisions + # self collisions for each species + explicit_fokker_planck_collisions_weak_form!(fvec_out.pdf,fvec_in.pdf,moments.charged.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.charged.dSdt, From f65699d7fe94911e719c3f428f9dda9a413d2d7f Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 15 Apr 2024 16:32:20 +0100 Subject: [PATCH 17/26] Make beam initial condition options more flexible. --- ...r-planck-relaxation-slowing-down-beam.toml | 5 +++++ moment_kinetics/src/initial_conditions.jl | 12 +++++------ moment_kinetics/src/input_structs.jl | 5 +++++ moment_kinetics/src/moment_kinetics_input.jl | 20 +++++++++++++------ 4 files changed, 30 insertions(+), 12 deletions(-) diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml index 7128f4173..849765df5 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-beam.toml @@ -27,6 +27,11 @@ 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 diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index a0849f17f..07d686c5e 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -691,9 +691,9 @@ function init_charged_pdf_over_density!(pdf, spec, composition, vpa, vperp, z, @. pdf[:,ivperp,iz] = (vpa.grid + 0.5*vpa.L)^spec.vpa_IC.monomial_degree end elseif spec.vpa_IC.initialization_option == "isotropic-beam" - v0 = 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2) # birth speed of beam - vslow = 0.2*v0 # spread of the beam in speed - v4norm = (vslow^2) * (v0^2) + 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 @@ -703,9 +703,9 @@ function init_charged_pdf_over_density!(pdf, spec, composition, vpa, vperp, z, @. pdf[:,:,iz] /= normfac end elseif spec.vpa_IC.initialization_option == "directed-beam" - vpa0 = 0.25*0.5*abs(vpa.L) # centre of beam in vpa - vperp0 = 0.5*abs(vperp.L) # centre of beam in vperp - vth0 = 0.05*sqrt(vperp.L^2 + (0.5*vpa.L)^2) # width of beam in v + 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 diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 0b25dcbce..7b6d496ea 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -209,6 +209,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 """ diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 65c61b8bc..06083c76c 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -449,19 +449,23 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) species.charged[is].z_IC.density_amplitude, species.charged[is].z_IC.density_phase, species.charged[is].z_IC.upar_amplitude, species.charged[is].z_IC.upar_phase, species.charged[is].z_IC.temperature_amplitude, species.charged[is].z_IC.temperature_phase, - species.charged[is].z_IC.monomial_degree) + species.charged[is].z_IC.monomial_degree, 0.0, 0.0, 0.0, 0.0) r_IC = initial_condition_input(species.charged[is].r_IC.initialization_option, species.charged[is].r_IC.width, species.charged[is].r_IC.wavenumber, species.charged[is].r_IC.density_amplitude, species.charged[is].r_IC.density_phase, species.charged[is].r_IC.upar_amplitude, species.charged[is].r_IC.upar_phase, species.charged[is].r_IC.temperature_amplitude, species.charged[is].r_IC.temperature_phase, - species.charged[is].r_IC.monomial_degree) + species.charged[is].r_IC.monomial_degree, 0.0, 0.0, 0.0, 0.0) vpa_IC = initial_condition_input(species.charged[is].vpa_IC.initialization_option, species.charged[is].vpa_IC.width, species.charged[is].vpa_IC.wavenumber, species.charged[is].vpa_IC.density_amplitude, species.charged[is].vpa_IC.density_phase, species.charged[is].vpa_IC.upar_amplitude, species.charged[is].vpa_IC.upar_phase, species.charged[is].vpa_IC.temperature_amplitude, - species.charged[is].vpa_IC.temperature_phase, species.charged[is].vpa_IC.monomial_degree) + species.charged[is].vpa_IC.temperature_phase, species.charged[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_charged_immutable[is] = species_parameters(species_type, species.charged[is].initial_temperature, species.charged[is].initial_density, z_IC, r_IC, vpa_IC) end @@ -473,19 +477,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 From 6e92bdef1d079b203867fd2e6bc9bdfe3c8ce4bb Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 17 Apr 2024 14:18:32 +0100 Subject: [PATCH 18/26] Additional plots of f(vpa,vperp) to aid diagnosis of slowing down simulations. --- .../src/plots_post_processing.jl | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) 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 54ba4e26c..e687ab1fa 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 @@ -2933,12 +2933,31 @@ function plot_charged_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r 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 From e7fc149d8feabae190532eee7ccf9469bcf9cb63 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 18 Apr 2024 14:27:32 +0100 Subject: [PATCH 19/26] Reinstate normalisation factor for alpha source. --- moment_kinetics/src/external_sources.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 02f94479d..c911ec047 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -530,7 +530,7 @@ function external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vp # Factor of 1/sqrt(π) (for 1V) or 1/π^(3/2) (for 2V/3V) is absorbed by the # normalisation of F pdf[ivpa,ivperp,iz,ir,is] += - this_prefactor * dummy_vpavperp[ivpa,ivperp] #/ normfac + this_prefactor * dummy_vpavperp[ivpa,ivperp] / normfac end end From bc63e2814c17891d663562481f87498a3a0068da Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 22 Apr 2024 11:29:45 +0100 Subject: [PATCH 20/26] Minor refactor of external sink term. --- moment_kinetics/src/external_sources.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index c911ec047..cdc47245b 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -546,14 +546,17 @@ function external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vp @loop_vperp_vpa ivperp ivpa begin v2 = vperp_grid[ivperp]^2 + vpa_grid[ivpa]^2 fac = 1.0/(sink_vth^2) - dummy_vpavperp[ivpa,ivperp] = (sink_strength/(sink_vth^3))*exp(-fac*v2) + dummy_vpavperp[ivpa,ivperp] = (1.0/(sink_vth^3))*exp(-fac*v2) end + # numerical correction to normalisation + normfac = get_density(dummy_vpavperp, vpa, vperp) + # println("sink norm", normfac) # add the source @loop_vperp_vpa ivperp ivpa begin # Factor of 1/sqrt(π) (for 1V) or 1/π^(3/2) (for 2V/3V) is absorbed by the # normalisation of F pdf[ivpa,ivperp,iz,ir,is] -= - dt * dummy_vpavperp[ivpa,ivperp] * pdf[ivpa,ivperp,iz,ir,is] + dt * sink_strength * dummy_vpavperp[ivpa,ivperp] * pdf[ivpa,ivperp,iz,ir,is] / normfac end end end From 8e5b91bcd5f67b52f762021b23053cf351a22581 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 22 Apr 2024 11:31:27 +0100 Subject: [PATCH 21/26] Correct the definition of the normalised collision frequency. Note that the code still only handles a single ion species so normalisation witth respect to charge number may have the be correctly stated later if variable charge numbers are added. --- moment_kinetics/src/fokker_planck.jl | 6 +++--- moment_kinetics/src/input_structs.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 9f9539f62..cae652911 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -81,8 +81,8 @@ nuii = 1.0 frequency_option = "manual" """ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) - # get reference collision frequency (note factor of 1/4 due to definition choices) - nuii_fkpl_default = 0.25*get_reference_collision_frequency(reference_params) + # get reference collision frequency (note factor of 1/2 due to definition choices) + nuii_fkpl_default = 0.5*get_reference_collision_frequency(reference_params) # read the input toml and specify a sensible default input_section = set_defaults_and_check_section!(toml_input, "fokker_planck_collisions", # begin default inputs (as kwargs) @@ -339,7 +339,7 @@ The result is stored in the array `fkpl_arrays.CC`. The normalised collision frequency is defined by ```math -\\tilde{\\nu}_{ss'} = \\frac{L_{\\mathrm{ref}}}{c_{\\mathrm{ref}}}\\frac{\\gamma_{ss'} n_\\mathrm{ref}}{2 m_s^2 c_\\mathrm{ref}^3} +\\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\$. diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 7b6d496ea..450b4c467 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -325,7 +325,7 @@ 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'} = (L/c_{ref}) * gamma_{ss'} n_{ref} / 2 (m_s)^2 (c_{ref})^3 + # nu_{ss'} = (L/c_{ref}) * 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 nuii::mk_float # option to determine if self collisions are used (for physics test) From 5944a81d71fa895f2b74a2026233c492e4a6dd9b Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 24 Apr 2024 09:27:49 +0100 Subject: [PATCH 22/26] Addition of "beam" and "beam-with-losses" options to external sources to model injection of particles around specific vpa0, vperp0, with a Gaussian shape function. In the "beam-with-losses" option, a sink is added to permit steady state solutions in a 0D2V simulation. --- moment_kinetics/src/external_sources.jl | 75 +++++++++++++++++++++++-- 1 file changed, 71 insertions(+), 4 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index cdc47245b..8d0365aab 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -48,15 +48,17 @@ function setup_external_sources!(input_dict, r, z) source_n=1.0, source_T=neutrals ? get(input_dict, "T_wall", 1.0) : 1.0, source_v0=0.0, # birth speed for "alphas" option - sink_strength=1.0, # strength of sink in "alphas-with-losses" option - sink_vth=0.0, # thermal speed for sink in "alphas-with-losses" option + source_vpa0=0.0, # birth vpa for "beam" option + source_vperp0=0.0, # birth vperp for "beam" option + sink_strength=1.0, # strength of sink in "alphas-with-losses" & "beam-with-losses" option + sink_vth=0.0, # thermal speed for sink in "alphas-with-losses" & "beam-with-losses" option r_profile="constant", r_width=1.0, r_relative_minimum=0.0, z_profile="constant", z_width=1.0, z_relative_minimum=0.0, - source_type="Maxwellian", # "energy", "alphas", "alphas-with-losses" + source_type="Maxwellian", # "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses" PI_density_controller_P=0.0, PI_density_controller_I=0.0, PI_density_target_amplitude=1.0, @@ -175,7 +177,7 @@ function setup_external_sources!(input_dict, r, z) PI_density_target_ir = nothing PI_density_target_iz = nothing PI_density_target_rank = nothing - elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses") + elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses") PI_density_target = nothing PI_controller_amplitude = nothing controller_source_profile = nothing @@ -564,6 +566,67 @@ function external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vp error("Unsupported combination in source_type=$(source_type) evolve_density=$(moments.evolve_density), " * "evolve_upar=$(moments.evolve_upar), evolve_ppar=$(moments.evolve_ppar)") end + elseif source_type == "beam" || source_type == "beam-with-losses" + begin_s_r_z_region() + source_vpa0 = ion_source_settings.source_vpa0 + source_vperp0 = ion_source_settings.source_vperp0 + if !(source_vpa0 > 1.0e-8) + error("source_vpa0=$source_vpa0 < 1.0e-8") + end + if !(source_vperp0 > 1.0e-8) + error("source_vperp0=$source_vperp0 < 1.0e-8") + end + dummy_vpavperp = scratch_dummy.dummy_vpavperp + if !moments.evolve_ppar && !moments.evolve_upar && !moments.evolve_density + @loop_s_r_z is ir iz begin + this_prefactor = dt * source_amplitude[iz,ir] + # first assign source to local scratch array + @loop_vperp_vpa ivperp ivpa begin + vth0 = sqrt(2.0*source_T) # sqrt(2 T / m), m = mref = 1 + v2 = ((vperp_grid[ivperp]-source_vperp0)^2 + (vpa_grid[ivpa]-source_vpa0)^2)/(vth0^2) + dummy_vpavperp[ivpa,ivperp] = (1.0/vth0^3)*exp(-v2) + end + # get the density for normalisation purposes + normfac = get_density(dummy_vpavperp, vpa, vperp) + # add the source + @loop_vperp_vpa ivperp ivpa begin + # Factor of 1/sqrt(π) (for 1V) or 1/π^(3/2) (for 2V/3V) is absorbed by the + # normalisation of F + pdf[ivpa,ivperp,iz,ir,is] += + this_prefactor * dummy_vpavperp[ivpa,ivperp] / normfac + end + end + + if source_type == "beam-with-losses" + sink_vth = ion_source_settings.sink_vth + sink_strength = ion_source_settings.sink_strength + if !(sink_vth > 1.0e-8) + error("sink_vth=$sink_vth < 1.0e-8") + end + # subtract a sink function representing the loss of slow ash particles + @loop_s_r_z is ir iz begin + # first assign sink to local scratch array + @loop_vperp_vpa ivperp ivpa begin + v2 = vperp_grid[ivperp]^2 + vpa_grid[ivpa]^2 + fac = 1.0/(sink_vth^2) + dummy_vpavperp[ivpa,ivperp] = (1.0/(sink_vth^3))*exp(-fac*v2) + end + # numerical correction to normalisation + normfac = get_density(dummy_vpavperp, vpa, vperp) + # println("sink norm", normfac) + # add the source + @loop_vperp_vpa ivperp ivpa begin + # Factor of 1/sqrt(π) (for 1V) or 1/π^(3/2) (for 2V/3V) is absorbed by the + # normalisation of F + pdf[ivpa,ivperp,iz,ir,is] -= + dt * sink_strength * dummy_vpavperp[ivpa,ivperp] * pdf[ivpa,ivperp,iz,ir,is] / normfac + end + end + end + else + error("Unsupported combination in source_type=$(source_type) evolve_density=$(moments.evolve_density), " + * "evolve_upar=$(moments.evolve_upar), evolve_ppar=$(moments.evolve_ppar)") + end else error("Unsupported source_type=$(source_type) ") end @@ -786,6 +849,10 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, # do nothing elseif ion_source_settings.source_type == "alphas-with-losses" # do nothing + elseif ion_source_settings.source_type == "beam" + # do nothing + elseif ion_source_settings.source_type == "beam-with-losses" + # do nothing else error("Unrecognised source_type=$(ion_source_settings.source_type)") end From 0b03e9bec6e934a2bd50b09ae9ee77063e3ae080 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 24 Apr 2024 16:33:53 +0100 Subject: [PATCH 23/26] Changes to make the Z of the evolved species an input parameter for the fokker-planck collision operator only. This adds Zi to the fokker planck input options. In the future, Zi could be inputted differently, but it will still be needed by the collision operator routines for calculating the correct collision frequency prefactor for each of the self and inter-species collision terms. --- moment_kinetics/src/fokker_planck.jl | 63 +++++++++++++++++++++++----- moment_kinetics/src/input_structs.jl | 15 +++++-- 2 files changed, 63 insertions(+), 15 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index cae652911..c7be07af8 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -91,11 +91,13 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) frequency_option = "manual", self_collisions = true, slowing_down_test = false, + #normalise_to_slowing_down_time = false, sd_density = 1.0, sd_temp = 0.01, - sd_q = 0.5, + sd_q = 1.0, sd_mi = 0.25, - sd_me = 0.25/1836.0) + sd_me = 0.25/1836.0, + Zi = 2.0) # ensure that the collision frequency is consistent with the input option frequency_option = input_section["frequency_option"] if frequency_option == "reference_parameters" @@ -114,6 +116,35 @@ 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 + #if input_section["normalise_to_slowing_down_time"] + # nuii = input_section["nuii"] + # input_section["nuii"] = nuii/nu_alphae + #end + #println(input_section["nuii"]) + end return fkpl_collisions_input(; input...) end @@ -234,14 +265,16 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ fkin = collisions.fkpl # masses charge numbers and collision frequencies - mref = 1.0 - Zref = 1.0 - nuref = fkin.nuii # generalise! + 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, fkin.sd_q] + Zsp = [fkin.sd_q, -1.0] # assume here that ne = sum_i n_i and that initial condition # for beam ions has unit density - densp = [fkin.sd_density, fkin.sd_density+1.0] + 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])] @@ -250,7 +283,7 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ @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,Zref,msp,Zsp,densp,uparsp,vthsp, + 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) @@ -296,7 +329,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 @@ -337,12 +372,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, @@ -513,7 +554,7 @@ function fokker_planck_collision_operator_weak_form_Maxwellian_Fsp!(ffs_in, dens = densp[isp] upar = uparsp[isp] vth = vthsp[isp] - ZZ = (Zsp[isp]/Zs)^2 # factor from gamma_ss' + 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) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 450b4c467..da4582467 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -324,27 +324,34 @@ 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'} = (L/c_{ref}) * 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 + #normalise_to_slowing_down_time::Bool + # Setting to force nuii -> nuii*(1/nu_alphae) so that tref = 1/nu_alphae 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 + # 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 """ From 939275a4526177d17017a988a935d7115eebd345 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 24 Apr 2024 17:19:47 +0100 Subject: [PATCH 24/26] Change default Zi for collisions to be unity. --- moment_kinetics/src/fokker_planck.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index a5759cd74..363c1e95e 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -97,7 +97,7 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) sd_q = 1.0, sd_mi = 0.25, sd_me = 0.25/1836.0, - Zi = 2.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" From a8d699bb520193fe4eb3d6e3dd673e8735d8c308 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 24 Apr 2024 17:20:54 +0100 Subject: [PATCH 25/26] Delete commented code. --- moment_kinetics/src/fokker_planck.jl | 6 ------ moment_kinetics/src/input_structs.jl | 2 -- 2 files changed, 8 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 363c1e95e..5b2922328 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -91,7 +91,6 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) frequency_option = "reference_parameters", self_collisions = true, slowing_down_test = false, - #normalise_to_slowing_down_time = false, sd_density = 1.0, sd_temp = 0.01, sd_q = 1.0, @@ -139,11 +138,6 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) println("vc3/cref^3 = $vc3") println("critical speed vc/cref = ",vc3^(1.0/3.0)) end - #if input_section["normalise_to_slowing_down_time"] - # nuii = input_section["nuii"] - # input_section["nuii"] = nuii/nu_alphae - #end - #println(input_section["nuii"]) end return fkpl_collisions_input(; input...) end diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index da4582467..6a19f1d1c 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -334,8 +334,6 @@ Base.@kwdef struct fkpl_collisions_input # 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 - #normalise_to_slowing_down_time::Bool - # Setting to force nuii -> nuii*(1/nu_alphae) so that tref = 1/nu_alphae frequency_option::String # "manual" # "reference_parameters" # options for fixed Maxwellian species in slowing down test operator # ion density - electron density determined from quasineutrality From cc91689eb597f84e2477b275de7cc085330b7abb Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 29 Apr 2024 17:23:56 +0100 Subject: [PATCH 26/26] Example slowing-down input file with source-sink behaviour, beam initial condition and physical parameters for alpha ion species. --- ...ation-slowing-down-alphas-source-sink.toml | 91 +++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml 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