diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index c5dbbcc32..f17bb33e4 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -18,7 +18,10 @@ jobs: - uses: julia-actions/cache@v1 - name: Install dependencies run: | - pip3 install --user matplotlib + # Version 3.9.0 of matplotlib causes an error with PyPlot.jl, so pin + # matplotlib version to 3.8.3 until this is fixed. See + # https://github.com/JuliaPy/PyPlot.jl/issues/582). + pip3 install --user "matplotlib==3.8.3" - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # Authenticate with GitHub Actions token @@ -33,6 +36,9 @@ jobs: version: '1.10' - name: Install dependencies run: | - pip3 install --user matplotlib + # Version 3.9.0 of matplotlib causes an error with PyPlot.jl, so pin + # matplotlib version to 3.8.3 until this is fixed. See + # https://github.com/JuliaPy/PyPlot.jl/issues/582). + pip3 install --user "matplotlib==3.8.3" - name: Build and deploy run: julia --project=docs/ docs/make-pdf.jl diff --git a/docs/src/zz_boundary_conditions.md b/docs/src/zz_boundary_conditions.md new file mode 100644 index 000000000..f10cbbc2a --- /dev/null +++ b/docs/src/zz_boundary_conditions.md @@ -0,0 +1,6 @@ +`boundary_conditions` +===================== + +```@autodocs +Modules = [moment_kinetics.boundary_conditions] +``` diff --git a/docs/src/zz_gyroaverages.md b/docs/src/zz_gyroaverages.md new file mode 100644 index 000000000..feb791dad --- /dev/null +++ b/docs/src/zz_gyroaverages.md @@ -0,0 +1,6 @@ +`gyroaverages` +============== + +```@autodocs +Modules = [moment_kinetics.gyroaverages] +``` diff --git a/docs/src/zz_lagrange_polynomials.md b/docs/src/zz_lagrange_polynomials.md new file mode 100644 index 000000000..af400823e --- /dev/null +++ b/docs/src/zz_lagrange_polynomials.md @@ -0,0 +1,6 @@ +`lagrange_polynomials` +====================== + +```@autodocs +Modules = [moment_kinetics.lagrange_polynomials] +``` diff --git a/moment_kinetics/src/calculus.jl b/moment_kinetics/src/calculus.jl index 31f07159b..fbff04a2f 100644 --- a/moment_kinetics/src/calculus.jl +++ b/moment_kinetics/src/calculus.jl @@ -14,6 +14,8 @@ using ..communication: block_rank using ..communication: _block_synchronize using ..looping +using LinearAlgebra + """ elementwise_derivative!(coord, f, adv_fac, spectral) elementwise_derivative!(coord, f, spectral) @@ -27,17 +29,6 @@ Result is stored in coord.scratch_2d. """ function elementwise_derivative! end -""" - elementwise_second_derivative!(coord, f, spectral) - -Generic function for element-by-element second derivatives. - -Note: no upwinding versions of second deriatives. - -Result is stored in coord.scratch_2d. -""" -function elementwise_second_derivative! end - """ derivative!(df, f, coord, adv_fac, spectral) @@ -73,7 +64,7 @@ function derivative!(df, f, coord, spectral::Union{null_spatial_dimension_info, return nothing end -function second_derivative!(d2f, f, coord, spectral) +function second_derivative!(d2f, f, coord, spectral; handle_periodic=true) # computes d^2f / d(coord)^2 # For spectral element methods, calculate second derivative by applying first # derivative twice, with special treatment for element boundaries @@ -152,40 +143,54 @@ is an input. """ function mass_matrix_solve! end -""" -Apply 'K-matrix' as part of a weak-form second derivative -""" -function elementwise_apply_Kmat! end - -function second_derivative!(d2f, f, coord, spectral::weak_discretization_info) +function second_derivative!(d2f, f, coord, spectral::weak_discretization_info; handle_periodic=true) # obtain the RHS of numerical weak-form of the equation # g = d^2 f / d coord^2, which is # M * g = K * f, with M the mass matrix and K an appropriate stiffness matrix # by multiplying by basis functions and integrating by parts - elementwise_apply_Kmat!(coord, f, spectral) - # map the RHS vector K * f from the elemental grid to the full grid; - # at element boundaries, use the average of K * f from neighboring elements. - derivative_elements_to_full_grid!(coord.scratch, coord.scratch_2d, coord) + mul!(coord.scratch, spectral.K_matrix, f) + + if handle_periodic && coord.bc == "periodic" + if coord.nrank > 1 + error("second_derivative!() cannot handle periodic boundaries for a " + * "distributed coordinate") + end + + coord.scratch[1] = 0.5 * (coord.scratch[1] + coord.scratch[end]) + coord.scratch[end] = coord.scratch[1] + end + # solve weak form matrix problem M * g = K * f to obtain g = d^2 f / d coord^2 + if coord.nrank > 1 + error("mass_matrix_solve!() does not support a " + * "distributed coordinate") + end mass_matrix_solve!(d2f, coord.scratch, spectral) end -""" -Apply 'L-matrix' as part of a weak-form Laplacian derivative -""" -function elementwise_apply_Lmat! end - function laplacian_derivative!(d2f, f, coord, spectral::weak_discretization_info) # for coord.name 'vperp' obtain the RHS of numerical weak-form of the equation # g = (1/coord) d/d coord ( coord d f / d coord ), which is # M * g = K * f, with M the mass matrix, and K an appropriate stiffness matrix, # by multiplying by basis functions and integrating by parts. # for all other coord.name, do exactly the same as second_derivative! above. - elementwise_apply_Lmat!(coord, f, spectral) - # map the RHS vector K * f from the elemental grid to the full grid; - # at element boundaries, use the average of K * f from neighboring elements. - derivative_elements_to_full_grid!(coord.scratch, coord.scratch_2d, coord) + mul!(coord.scratch, spectral.L_matrix, f) + + if handle_periodic && coord.bc == "periodic" + if coord.nrank > 1 + error("second_derivative!() cannot handle periodic boundaries for a " + * "distributed coordinate") + end + + coord.scratch[1] = 0.5 * (coord.scratch[1] + coord.scratch[end]) + coord.scratch[end] = coord.scratch[1] + end + # solve weak form matrix problem M * g = K * f to obtain g = d^2 f / d coord^2 + if coord.nrank > 1 + error("mass_matrix_solve!() does not support a " + * "distributed coordinate") + end mass_matrix_solve!(d2f, coord.scratch, spectral) end diff --git a/moment_kinetics/src/chebyshev.jl b/moment_kinetics/src/chebyshev.jl index e3c789aec..ad8139120 100644 --- a/moment_kinetics/src/chebyshev.jl +++ b/moment_kinetics/src/chebyshev.jl @@ -18,7 +18,7 @@ using ..array_allocation: allocate_float, allocate_complex using ..clenshaw_curtis: clenshawcurtisweights import ..calculus: elementwise_derivative! using ..communication -import ..interpolation: interpolate_to_grid_1d! +import ..interpolation: single_element_interpolate! using ..moment_kinetics_structs: discretization_info """ @@ -465,88 +465,8 @@ function chebyshev_spectral_derivative!(df,f) end end -""" -Interpolation from a regular grid to a 1d grid with arbitrary spacing - -Arguments ---------- -result : Array{mk_float, 1} - Array to be overwritten with the result of the interpolation -new_grid : Array{mk_float, 1} - Grid of points to interpolate `coord` to -f : Array{mk_float} - Field to be interpolated -coord : coordinate - `coordinate` struct giving the coordinate along which f varies -chebyshev : chebyshev_info - struct containing information for Chebyshev transforms - -Note that this routine does not support Gauss-Chebyshev-Radau elements -""" -function interpolate_to_grid_1d!(result, newgrid, f, coord, chebyshev::chebyshev_info) - # define local variable nelement for convenience - nelement = coord.nelement_local - # check array bounds - @boundscheck nelement == size(chebyshev.lobatto.f,2) || throw(BoundsError(chebyshev.lobatto.f)) - - n_new = size(newgrid)[1] - # Find which points belong to which element. - # kstart[j] contains the index of the first point in newgrid that is within element - # j, and kstart[nelement+1] is n_new if the last point is within coord.grid, or the - # index of the first element outside coord.grid otherwise. - # Assumes points in newgrid are sorted. - # May not be the most efficient algorithm. - # Find start/end points for each element, storing their indices in kstart - kstart = Vector{mk_int}(undef, nelement+1) - # set the starting index by finding the start of coord.grid - kstart[1] = searchsortedfirst(newgrid, coord.grid[1]) - # check to see if any of the newgrid points are to the left of the first grid point - for j ∈ 1:kstart[1]-1 - # if the new grid location is outside the bounds of the original grid, - # extrapolate f with Gaussian-like decay beyond the domain - result[j] = f[1] * exp(-(coord.grid[1] - newgrid[j])^2) - end - @inbounds for j ∈ 1:nelement - # Search from kstart[j] to try to speed up the sort, but means result of - # searchsortedfirst() is offset by kstart[j]-1 from the beginning of newgrid. - kstart[j+1] = kstart[j] - 1 + @views searchsortedfirst(newgrid[kstart[j]:end], coord.grid[coord.imax[j]]) - end - - # First element includes both boundary points, while all others have only one (to - # avoid duplication), so calculate the first element outside the loop. - if kstart[1] < kstart[2] - imin = coord.imin[1] - imax = coord.imax[1] - kmin = kstart[1] - kmax = kstart[2] - 1 - @views chebyshev_interpolate_single_element!(result[kmin:kmax], - newgrid[kmin:kmax], - f[imin:imax], - imin, imax, coord, chebyshev.lobatto) - end - @inbounds for j ∈ 2:nelement - kmin = kstart[j] - kmax = kstart[j+1] - 1 - if kmin <= kmax - imin = coord.imin[j] - 1 - imax = coord.imax[j] - @views chebyshev_interpolate_single_element!(result[kmin:kmax], - newgrid[kmin:kmax], - f[imin:imax], - imin, imax, coord, chebyshev.lobatto) - end - end - - for k ∈ kstart[nelement+1]:n_new - result[k] = f[end] * exp(-(newgrid[k] - coord.grid[end])^2) - end - - return nothing -end - -""" -""" -function chebyshev_interpolate_single_element!(result, newgrid, f, imin, imax, coord, chebyshev::chebyshev_base_info) +function single_element_interpolate!(result, newgrid, f, imin, imax, ielement, coord, + chebyshev::chebyshev_base_info) # Temporary buffer to store Chebyshev coefficients cheby_f = chebyshev.df diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl index 95ce99da6..5d5531d84 100644 --- a/moment_kinetics/src/coordinates.jl +++ b/moment_kinetics/src/coordinates.jl @@ -107,6 +107,12 @@ struct coordinate{T <: AbstractVector{mk_float}} element_spacing_option::String # list of element boundaries element_boundaries::Array{mk_float,1} + # Does the coordinate use a 'Radau' discretization for the first element? + radau_first_element::Bool + # 'Other' nodes where the j'th Lagrange polynomial (which is 1 at x[j]) is equal to 0 + other_nodes::Array{mk_float,3} + # One over the denominators of the Lagrange polynomials + one_over_denominator::Array{mk_float,2} end """ @@ -115,7 +121,7 @@ setup the coordinate grid, and populate the coordinate structure containing all of this information """ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing, - ignore_MPI=false, init_YY::Bool=true) + ignore_MPI=false, collision_operator_dim::Bool=true) # total number of grid points is ngrid for the first element # plus ngrid-1 unique points for each additional element due # to the repetition of a point at the element boundary @@ -136,8 +142,9 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing element_scale, element_shift = set_element_scale_and_shift(input.nelement_global, input.nelement_local, input.irank, element_boundaries) # initialize the grid and the integration weights associated with the grid # also obtain the Chebyshev theta grid and spacing if chosen as discretization option - grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_local, n_global, n_local, input.irank, input.L, element_scale, element_shift, - imin, imax, igrid, input.discretization, input.name) + grid, wgts, uniform_grid, radau_first_element = init_grid(input.ngrid, + input.nelement_local, n_global, n_local, input.irank, input.L, element_scale, + element_shift, imin, imax, igrid, input.discretization, input.name) # calculate the widths of the cells between neighboring grid points cell_width = grid_spacing(grid, n_local) # duniform_dgrid is the local derivative of the uniform grid with respect to @@ -187,13 +194,37 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing local_io_range = 1 : n_local-1 global_io_range = input.irank*(n_local-1)+1 : (input.irank+1)*(n_local-1) end + + # Precompute some values for Lagrange polynomial evaluation + other_nodes = allocate_float(input.ngrid-1, input.ngrid, input.nelement_local) + one_over_denominator = allocate_float(input.ngrid, input.nelement_local) + for ielement ∈ 1:input.nelement_local + if ielement == 1 + this_imin = imin[ielement] + else + this_imin = imin[ielement] - 1 + end + this_imax = imax[ielement] + this_grid = grid[this_imin:this_imax] + for j ∈ 1:input.ngrid + @views other_nodes[1:j-1,j,ielement] .= this_grid[1:j-1] + @views other_nodes[j:end,j,ielement] .= this_grid[j+1:end] + + if input.ngrid == 1 + one_over_denominator[j,ielement] = 1.0 + else + one_over_denominator[j,ielement] = 1.0 / prod(this_grid[j] - n for n ∈ @view other_nodes[:,j,ielement]) + end + end + end + coord = coordinate(input.name, n_global, n_local, input.ngrid, input.nelement_global, input.nelement_local, input.nrank, input.irank, input.L, grid, 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_shared, scratch_shared2, 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, - element_boundaries) + element_boundaries, radau_first_element, other_nodes, one_over_denominator) if coord.n == 1 && occursin("v", coord.name) spectral = null_velocity_dimension_info() @@ -211,7 +242,7 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing elseif input.discretization == "gausslegendre_pseudospectral" # create arrays needed for explicit GaussLegendre pseudospectral treatment in this # coordinate and create the matrices for differentiation - spectral = setup_gausslegendre_pseudospectral(coord,init_YY=init_YY) + spectral = setup_gausslegendre_pseudospectral(coord, collision_operator_dim=collision_operator_dim) # obtain the local derivatives of the uniform grid with respect to the used grid derivative!(coord.duniform_dgrid, coord.uniform_grid, coord, spectral) else @@ -283,6 +314,7 @@ function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_s imin, imax, igrid, discretization, name) uniform_grid = equally_spaced_grid(n_global, n_local, irank, L) uniform_grid_shifted = equally_spaced_grid_shifted(n_global, n_local, irank, L) + radau_first_element = false if n_global == 1 grid = allocate_float(n_local) grid[1] = 0.0 @@ -299,6 +331,7 @@ function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_s grid, wgts = scaled_chebyshev_radau_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax, irank) wgts = 2.0 .* wgts .* grid # to include 2 vperp in jacobian of integral # see note above on normalisation + radau_first_element = true else # initialize chebyshev grid defined on [-L/2,L/2] # with n grid points chosen to facilitate @@ -314,6 +347,7 @@ function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_s grid, wgts = scaled_gauss_legendre_radau_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax, irank) wgts = 2.0 .* wgts .* grid # to include 2 vperp in jacobian of integral # see note above on normalisation + radau_first_element = true else grid, wgts = scaled_gauss_legendre_lobatto_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax) end @@ -336,7 +370,7 @@ function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_s error("discretization option '$discretization' unrecognized") end # return the locations of the grid points - return grid, wgts, uniform_grid + return grid, wgts, uniform_grid, radau_first_element end """ diff --git a/moment_kinetics/src/derivatives.jl b/moment_kinetics/src/derivatives.jl index 409274b0f..e85e91158 100644 --- a/moment_kinetics/src/derivatives.jl +++ b/moment_kinetics/src/derivatives.jl @@ -11,7 +11,7 @@ module derivatives export derivative_r!, derivative_r_chrg!, derivative_r_ntrl! export derivative_z!, derivative_z_chrg!, derivative_z_ntrl! -using ..calculus: derivative!, reconcile_element_boundaries_MPI! +using ..calculus: derivative!, second_derivative!, reconcile_element_boundaries_MPI! using ..type_definitions: mk_float using ..looping @@ -266,6 +266,251 @@ function derivative_z!(dfdz::AbstractArray{mk_float,6}, f::AbstractArray{mk_floa end +""" +Centered derivatives +df2/dr2 group of rountines for +fields & moments -> [z,r] +dfns (ion) -> [vpa,vperp,z,r,s] +dfns (neutrals) -> [vz,vr,vzeta,z,r,sn] +""" + +#d2f/dr2 +#2D version for f[z,r] -> Er, Ez, phi +function second_derivative_r!(d2fdr2::AbstractArray{mk_float,2}, f::AbstractArray{mk_float,2}, + d2fdr2_lower_endpoints::AbstractArray{mk_float,1}, + d2fdr2_upper_endpoints::AbstractArray{mk_float,1}, + r_receive_buffer1::AbstractArray{mk_float,1}, + r_receive_buffer2::AbstractArray{mk_float,1}, r_spectral, r) + + begin_z_region() + + # differentiate f w.r.t r + @loop_z iz begin + @views second_derivative!(d2fdr2[iz,:], f[iz,:], r, r_spectral) + # get external endpoints to reconcile via MPI + d2fdr2_lower_endpoints[iz] = r.scratch_2d[1,1] + d2fdr2_upper_endpoints[iz] = r.scratch_2d[end,end] + end + # now reconcile element boundaries across + # processes with large message involving all other dimensions + if r.nelement_local < r.nelement_global + reconcile_element_boundaries_MPI!(d2fdr2, d2fdr2_lower_endpoints, + d2fdr2_upper_endpoints, r_receive_buffer1, + r_receive_buffer2, r) + end +end + +#d2f/dr2 +#3D version for f[s,z,r] -> moments n, u, T etc +function second_derivative_r!(d2fdr2::AbstractArray{mk_float,3}, f::AbstractArray{mk_float,3}, + d2fdr2_lower_endpoints::AbstractArray{mk_float,2}, + d2fdr2_upper_endpoints::AbstractArray{mk_float,2}, + r_receive_buffer1::AbstractArray{mk_float,2}, + r_receive_buffer2::AbstractArray{mk_float,2}, r_spectral, r; neutrals=false) + + # differentiate f w.r.t r + if neutrals + @loop_sn_z isn iz begin + @views second_derivative!(d2fdr2[iz,:,isn], f[iz,:,isn], r, r_spectral) + # get external endpoints to reconcile via MPI + d2fdr2_lower_endpoints[iz,isn] = r.scratch_2d[1,1] + d2fdr2_upper_endpoints[iz,isn] = r.scratch_2d[end,end] + end + else + @loop_s_z is iz begin + @views second_derivative!(d2fdr2[iz,:,is], f[iz,:,is], r, r_spectral) + # get external endpoints to reconcile via MPI + d2fdr2_lower_endpoints[iz,is] = r.scratch_2d[1,1] + d2fdr2_upper_endpoints[iz,is] = r.scratch_2d[end,end] + end + end + + # Sometimes an array might contain no data (e.g. if n_neutral_species=0). Then don't + # need to reconcile boundaries + if length(d2fdr2) > 0 + # now reconcile element boundaries across + # processes with large message involving all other dimensions + if r.nelement_local < r.nelement_global + reconcile_element_boundaries_MPI!(d2fdr2, d2fdr2_lower_endpoints, + d2fdr2_upper_endpoints, r_receive_buffer1, + r_receive_buffer2, r) + end + end +end + +#d2f/dr2 +#5D version for f[vpa,vperp,z,r,s] -> ion particle dfn +function second_derivative_r!(d2fdr2::AbstractArray{mk_float,5}, f::AbstractArray{mk_float,5}, + d2fdr2_lower_endpoints::AbstractArray{mk_float,4}, + d2fdr2_upper_endpoints::AbstractArray{mk_float,4}, + r_receive_buffer1::AbstractArray{mk_float,4}, + r_receive_buffer2::AbstractArray{mk_float,4}, r_spectral, r) + + begin_s_z_vperp_vpa_region() + + # differentiate f w.r.t r + @loop_s_z_vperp_vpa is iz ivperp ivpa begin + @views second_derivative!(d2fdr2[ivpa,ivperp,iz,:,is], f[ivpa,ivperp,iz,:,is], r, r_spectral) + # get external endpoints to reconcile via MPI + d2fdr2_lower_endpoints[ivpa,ivperp,iz,is] = r.scratch_2d[1,1] + d2fdr2_upper_endpoints[ivpa,ivperp,iz,is] = r.scratch_2d[end,end] + end + # now reconcile element boundaries across + # processes with large message involving all other dimensions + if r.nelement_local < r.nelement_global + reconcile_element_boundaries_MPI!(d2fdr2, d2fdr2_lower_endpoints, + d2fdr2_upper_endpoints, r_receive_buffer1, + r_receive_buffer2, r) + end +end + +#6D version for f[vz,vz,vzeta,z,r,sn] -> neutral particle dfn (species indexing taken outside this loop) +function second_derivative_r!(d2fdr2::AbstractArray{mk_float,6}, f::AbstractArray{mk_float,6}, + d2fdr2_lower_endpoints::AbstractArray{mk_float,5}, + d2fdr2_upper_endpoints::AbstractArray{mk_float,5}, + r_receive_buffer1::AbstractArray{mk_float,5}, + r_receive_buffer2::AbstractArray{mk_float,5}, r_spectral, r) + + begin_sn_z_vzeta_vr_vz_region() + + # differentiate f w.r.t r + @loop_sn_z_vzeta_vr_vz isn iz ivzeta ivr ivz begin + @views second_derivative!(d2fdr2[ivz,ivr,ivzeta,iz,:,isn], f[ivz,ivr,ivzeta,iz,:,isn], r, r_spectral) + # get external endpoints to reconcile via MPI + d2fdr2_lower_endpoints[ivz,ivr,ivzeta,iz,isn] = r.scratch_2d[1,1] + d2fdr2_upper_endpoints[ivz,ivr,ivzeta,iz,isn] = r.scratch_2d[end,end] + end + # now reconcile element boundaries across + # processes with large message involving all other dimensions + if r.nelement_local < r.nelement_global + reconcile_element_boundaries_MPI!(d2fdr2, d2fdr2_lower_endpoints, + d2fdr2_upper_endpoints, r_receive_buffer1, + r_receive_buffer2, r) + end +end + +""" +Centered derivatives +df/dz group of rountines for +fields & moments -> [z,r] +dfns (ion) -> [vpa,vperp,z,r,s] +dfns (neutrals) -> [vz,vr,vzeta,z,r,sn] +""" + +#d2f/dz2 +#2D version for f[z,r] -> Er, Ez, phi +function second_derivative_z!(d2fdz2::AbstractArray{mk_float,2}, f::AbstractArray{mk_float,2}, + d2fdz2_lower_endpoints::AbstractArray{mk_float,1}, + d2fdz2_upper_endpoints::AbstractArray{mk_float,1}, + z_send_buffer::AbstractArray{mk_float,1}, + z_receive_buffer::AbstractArray{mk_float,1}, z_spectral, z) + + begin_r_region() + + # differentiate f w.r.t z + @loop_r ir begin + @views second_derivative!(d2fdz2[:,ir], f[:,ir], z, z_spectral) + # get external endpoints to reconcile via MPI + d2fdz2_lower_endpoints[ir] = z.scratch_2d[1,1] + d2fdz2_upper_endpoints[ir] = z.scratch_2d[end,end] + end + # now reconcile element boundaries across + # processes with large message involving all y + if z.nelement_local < z.nelement_global + reconcile_element_boundaries_MPI!(d2fdz2, d2fdz2_lower_endpoints, + d2fdz2_upper_endpoints, z_send_buffer, + z_receive_buffer, z) + end +end + +#d2f/dz2 +#3D version for f[z,r] -> moments n, u, T etc +function second_derivative_z!(d2fdz2::AbstractArray{mk_float,3}, f::AbstractArray{mk_float,3}, + d2fdz2_lower_endpoints::AbstractArray{mk_float,2}, + d2fdz2_upper_endpoints::AbstractArray{mk_float,2}, + z_send_buffer::AbstractArray{mk_float,2}, + z_receive_buffer::AbstractArray{mk_float,2}, z_spectral, z; neutrals=false) + + # differentiate f w.r.t z + if neutrals + @loop_sn_r isn ir begin + @views second_derivative!(d2fdz2[:,ir,isn], f[:,ir,isn], z, z_spectral) + # get external endpoints to reconcile via MPI + d2fdz2_lower_endpoints[ir,isn] = z.scratch_2d[1,1] + d2fdz2_upper_endpoints[ir,isn] = z.scratch_2d[end,end] + end + else + @loop_s_r is ir begin + @views second_derivative!(d2fdz2[:,ir,is], f[:,ir,is], z, z_spectral) + # get external endpoints to reconcile via MPI + d2fdz2_lower_endpoints[ir,is] = z.scratch_2d[1,1] + d2fdz2_upper_endpoints[ir,is] = z.scratch_2d[end,end] + end + end + + # Sometimes an array might contain no data (e.g. if n_neutral_species=0). Then don't + # need to reconcile boundaries + if length(d2fdz2) > 0 + # now reconcile element boundaries across + # processes with large message involving all y + if z.nelement_local < z.nelement_global + reconcile_element_boundaries_MPI!(d2fdz2, d2fdz2_lower_endpoints, + d2fdz2_upper_endpoints, z_send_buffer, + z_receive_buffer, z) + end + end +end + +#5D version for f[vpa,vperp,z,r,s] -> dfn ions +function second_derivative_z!(d2fdz2::AbstractArray{mk_float,5}, f::AbstractArray{mk_float,5}, + d2fdz2_lower_endpoints::AbstractArray{mk_float,4}, + d2fdz2_upper_endpoints::AbstractArray{mk_float,4}, + z_send_buffer::AbstractArray{mk_float,4}, + z_receive_buffer::AbstractArray{mk_float,4}, z_spectral, z) + + begin_s_r_vperp_vpa_region() + + # differentiate f w.r.t z + @loop_s_r_vperp_vpa is ir ivperp ivpa begin + @views second_derivative!(d2fdz2[ivpa,ivperp,:,ir,is], f[ivpa,ivperp,:,ir,is], z, z_spectral) + # get external endpoints to reconcile via MPI + d2fdz2_lower_endpoints[ivpa,ivperp,ir,is] = z.scratch_2d[1,1] + d2fdz2_upper_endpoints[ivpa,ivperp,ir,is] = z.scratch_2d[end,end] + end + # now reconcile element boundaries across + # processes with large message involving all y + if z.nelement_local < z.nelement_global + reconcile_element_boundaries_MPI!(d2fdz2, d2fdz2_lower_endpoints, + d2fdz2_upper_endpoints, z_send_buffer, + z_receive_buffer, z) + end +end + +#6D version for f[vz,vr,vzeta,z,r] -> dfn neutral particles +function second_derivative_z!(d2fdz2::AbstractArray{mk_float,6}, f::AbstractArray{mk_float,6}, + d2fdz2_lower_endpoints::AbstractArray{mk_float,5}, + d2fdz2_upper_endpoints::AbstractArray{mk_float,5}, + z_send_buffer::AbstractArray{mk_float,5}, + z_receive_buffer::AbstractArray{mk_float,5}, z_spectral, z) + + begin_sn_r_vzeta_vr_vz_region() + + # differentiate f w.r.t z + @loop_sn_r_vzeta_vr_vz isn ir ivzeta ivr ivz begin + @views second_derivative!(d2fdz2[ivz,ivr,ivzeta,:,ir,isn], f[ivz,ivr,ivzeta,:,ir,isn], z, z_spectral) + # get external endpoints to reconcile via MPI + d2fdz2_lower_endpoints[ivz,ivr,ivzeta,ir,isn] = z.scratch_2d[1,1] + d2fdz2_upper_endpoints[ivz,ivr,ivzeta,ir,isn] = z.scratch_2d[end,end] + end + # now reconcile element boundaries across + # processes with large message involving all y + if z.nelement_local < z.nelement_global + reconcile_element_boundaries_MPI!(d2fdz2, d2fdz2_lower_endpoints, + d2fdz2_upper_endpoints, z_send_buffer, + z_receive_buffer, z) + end +end + """ Upwind derivatives df/dr group of rountines for diff --git a/moment_kinetics/src/finite_differences.jl b/moment_kinetics/src/finite_differences.jl index 17f3c58c2..ae18c7d7d 100644 --- a/moment_kinetics/src/finite_differences.jl +++ b/moment_kinetics/src/finite_differences.jl @@ -3,8 +3,8 @@ module finite_differences using ..type_definitions: mk_float -import ..calculus: elementwise_derivative!, elementwise_second_derivative!, - second_derivative! +import ..calculus: elementwise_derivative!, second_derivative!, + derivative_elements_to_full_grid! using ..moment_kinetics_structs: discretization_info """ diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index ca8697628..300e0ca6e 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -39,6 +39,7 @@ using ..array_allocation: allocate_float, allocate_shared_float using ..calculus: derivative! using ..communication using ..communication: MPISharedArray, global_rank +using ..lagrange_polynomials: lagrange_poly, lagrange_poly_optimised using ..looping using moment_kinetics.gauss_legendre: get_QQ_local! using Dates @@ -506,28 +507,6 @@ function get_nodes(coord,iel) nodes = coord.grid[imin:imax] return nodes end -""" -Lagrange polynomial -args: -j - index of l_j from list of nodes -x_nodes - array of x node values -x - point where interpolated value is returned -""" -function lagrange_poly(j,x_nodes,x) - # get number of nodes - n = size(x_nodes,1) - # location where l(x0) = 1 - x0 = x_nodes[j] - # evaluate polynomial - poly = 1.0 - for i in 1:j-1 - poly *= (x - x_nodes[i])/(x0 - x_nodes[i]) - end - for i in j+1:n - poly *= (x - x_nodes[i])/(x0 - x_nodes[i]) - end - return poly -end # Function to get the local integration grid and quadrature weights # to integrate a 1D element in the 2D representation of the @@ -667,12 +646,16 @@ end # `ellipe(k) = \int^{\pi/2}\_0 \frac{1}{\sqrt{ 1 - m \sin^2(\theta)}} d \theta` function local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, - nquad_vpa,ielement_vpa,vpa_nodes,vpa, # info about primed vpa grids - nquad_vperp,ielement_vperp,vperp_nodes,vperp, # info about primed vperp grids + nquad_vpa,ielement_vpa,vpa, # info about primed vpa grids + nquad_vperp,ielement_vperp,vperp, # info about primed vperp grids x_vpa, w_vpa, x_vperp, w_vperp, # points and weights for primed (source) grids vpa_val, vperp_val) # values and indices for unprimed (field) grids for igrid_vperp in 1:vperp.ngrid + vperp_other_nodes = @view vperp.other_nodes[:,igrid_vperp,ielement_vperp] + vperp_one_over_denominator = vperp.one_over_denominator[igrid_vperp,ielement_vperp] for igrid_vpa in 1:vpa.ngrid + vpa_other_nodes = @view vpa.other_nodes[:,igrid_vpa,ielement_vpa] + vpa_one_over_denominator = vpa.one_over_denominator[igrid_vpa,ielement_vpa] # get grid index for point on full grid ivpap = vpa.igrid_full[igrid_vpa,ielement_vpa] ivperpp = vperp.igrid_full[igrid_vperp,ielement_vperp] @@ -700,8 +683,12 @@ function local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights, H_elliptic_integral_factor = 2.0*ellipk_mm/(pi*prefac) H1_elliptic_integral_factor = -(2.0/(pi*prefac))*( (mm-2.0)*(ellipk_mm/mm) + (2.0*ellipe_mm/mm) ) H2_elliptic_integral_factor = (2.0/(pi*prefac))*( (3.0*mm^2 - 8.0*mm + 8.0)*(ellipk_mm/(3.0*mm^2)) + (4.0*mm - 8.0)*ellipe_mm/(3.0*mm^2) ) - lagrange_poly_vpa = lagrange_poly(igrid_vpa,vpa_nodes,x_kvpa) - lagrange_poly_vperp = lagrange_poly(igrid_vperp,vperp_nodes,x_kvperp) + lagrange_poly_vpa = lagrange_poly_optimised(vpa_other_nodes, + vpa_one_over_denominator, + x_kvpa) + lagrange_poly_vperp = lagrange_poly_optimised(vperp_other_nodes, + vperp_one_over_denominator, + x_kvperp) (G0_weights[ivpap,ivperpp] += lagrange_poly_vpa*lagrange_poly_vperp* @@ -762,8 +749,8 @@ function loop_over_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weights,H2_ vpa_min, vpa_max = vpa_nodes[1], vpa_nodes[end] nquad_vpa = get_scaled_x_w_no_divergences!(x_vpa, w_vpa, x_legendre, w_legendre, vpa_min, vpa_max) local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, - nquad_vpa,ielement_vpap,vpa_nodes,vpa, - nquad_vperp,ielement_vperpp,vperp_nodes,vperp, + nquad_vpa,ielement_vpap,vpa, + nquad_vperp,ielement_vperpp,vperp, x_vpa, w_vpa, x_vperp, w_vperp, vpa_val, vperp_val) end @@ -776,8 +763,8 @@ function loop_over_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weights,H2_ #nquad_vpa = get_scaled_x_w_no_divergences!(x_vpa, w_vpa, x_legendre, w_legendre, vpa_min, vpa_max) nquad_vpa = get_scaled_x_w_with_divergences!(x_vpa, w_vpa, x_legendre, w_legendre, x_laguerre, w_laguerre, vpa_min, vpa_max, vpa_nodes, igrid_vpa, vpa_val) local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, - nquad_vpa,ielement_vpap,vpa_nodes,vpa, - nquad_vperp,ielement_vperpp,vperp_nodes,vperp, + nquad_vpa,ielement_vpap,vpa, + nquad_vperp,ielement_vperpp,vperp, x_vpa, w_vpa, x_vperp, w_vperp, vpa_val, vperp_val) end @@ -788,8 +775,8 @@ function loop_over_vpa_elements!(G0_weights,G1_weights,H0_weights,H1_weights,H2_ vpa_min, vpa_max = vpa_nodes[1], vpa_nodes[end] nquad_vpa = get_scaled_x_w_no_divergences!(x_vpa, w_vpa, x_legendre, w_legendre, vpa_min, vpa_max) local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, - nquad_vpa,ielement_vpap,vpa_nodes,vpa, - nquad_vperp,ielement_vperpp,vperp_nodes,vperp, + nquad_vpa,ielement_vpap,vpa, + nquad_vperp,ielement_vperpp,vperp, x_vpa, w_vpa, x_vperp, w_vperp, vpa_val, vperp_val) @@ -809,8 +796,8 @@ function loop_over_vpa_elements_no_divergences!(G0_weights,G1_weights,H0_weights vpa_min, vpa_max = vpa_nodes[1], vpa_nodes[end] nquad_vpa = get_scaled_x_w_no_divergences!(x_vpa, w_vpa, x_legendre, w_legendre, vpa_min, vpa_max) local_element_integration!(G0_weights,G1_weights,H0_weights,H1_weights,H2_weights,H3_weights, - nquad_vpa,ielement_vpap,vpa_nodes,vpa, - nquad_vperp,ielement_vperpp,vperp_nodes,vperp, + nquad_vpa,ielement_vpap,vpa, + nquad_vperp,ielement_vperpp,vperp, x_vpa, w_vpa, x_vperp, w_vperp, vpa_val, vperp_val) @@ -2387,6 +2374,25 @@ function interpolate_2D_vspace!(pdf_out,pdf_in,vpa,vperp,scalefac) end return nothing end +# Alternative version that should be faster - to be tested +#function interpolate_2D_vspace!(pdf_out, pdf_in, vpa, vpa_spectral, vperp, vperp_spectral, +# scalefac, pdf_buffer) +# newgrid_vperp = vperp.scratch .= scalefac .* vperp.grid +# newgrid_vpa = vpa.scratch .= scalefac .* vpa.grid +# +# begin_anyv_vpa_region() +# @loop_vpa ivpa begin +# @views interpolate_to_grid_1d!(pdf_buffer[ivpa,:], newgrid_vperp, +# pdf_in[ivpa,:], vperp, vperp_spectral) +# end +# +# begin_anyv_vperp_region() +# @loop_vperp ivperp begin +# @views interpolate_to_grid_1d!(pdf_out[:,ivperp], newgrid_vpa, +# pdf_buffer[:,ivperp], vpa, vpa_spectral) + +# end +#end """ function to find the element in which x sits diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl index c93d263c2..24a0b925f 100644 --- a/moment_kinetics/src/gauss_legendre.jl +++ b/moment_kinetics/src/gauss_legendre.jl @@ -27,8 +27,9 @@ using LinearAlgebra: mul!, lu, LU using SparseArrays: sparse, AbstractSparseArray using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float -import ..calculus: elementwise_derivative!, elementwise_apply_Kmat!, - elementwise_apply_Lmat!, mass_matrix_solve! +import ..calculus: elementwise_derivative!, mass_matrix_solve! +import ..interpolation: single_element_interpolate! +using ..lagrange_polynomials: lagrange_poly_optimised using ..moment_kinetics_structs: weak_discretization_info @@ -81,42 +82,48 @@ struct gausslegendre_base_info Y31::Array{mk_float,3} end -struct gausslegendre_info <: weak_discretization_info +struct gausslegendre_info{TSparse, TLU} <: weak_discretization_info lobatto::gausslegendre_base_info radau::gausslegendre_base_info # global (1D) mass matrix mass_matrix::Array{mk_float,2} # global (1D) weak derivative matrix #S_matrix::Array{mk_float,2} - S_matrix::AbstractSparseArray{mk_float,Ti,2} where Ti + S_matrix::TSparse # global (1D) weak second derivative matrix - K_matrix::Array{mk_float,2} + K_matrix::TSparse # global (1D) weak Laplacian derivative matrix - L_matrix::Array{mk_float,2} + L_matrix::TSparse # global (1D) LU object - mass_matrix_lu::T where T + mass_matrix_lu::TLU # dummy matrix for local operators Qmat::Array{mk_float,2} end -function setup_gausslegendre_pseudospectral(coord;init_YY=true) - lobatto = setup_gausslegendre_pseudospectral_lobatto(coord,init_YY=init_YY) - radau = setup_gausslegendre_pseudospectral_radau(coord,init_YY=init_YY) +function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true) + lobatto = setup_gausslegendre_pseudospectral_lobatto(coord,collision_operator_dim=collision_operator_dim) + radau = setup_gausslegendre_pseudospectral_radau(coord,collision_operator_dim=collision_operator_dim) + + if collision_operator_dim + S_matrix = allocate_float(coord.n,coord.n) + setup_global_weak_form_matrix!(S_matrix, lobatto, radau, coord, "S") + else + S_matrix = allocate_float(0, 0) + end mass_matrix = allocate_float(coord.n,coord.n) - S_matrix = allocate_float(coord.n,coord.n) K_matrix = allocate_float(coord.n,coord.n) L_matrix = allocate_float(coord.n,coord.n) - + setup_global_weak_form_matrix!(mass_matrix, lobatto, radau, coord, "M") - setup_global_weak_form_matrix!(S_matrix, lobatto, radau, coord, "S") setup_global_weak_form_matrix!(K_matrix, lobatto, radau, coord, "K_with_BC_terms") setup_global_weak_form_matrix!(L_matrix, lobatto, radau, coord, "L_with_BC_terms") mass_matrix_lu = lu(sparse(mass_matrix)) Qmat = allocate_float(coord.ngrid,coord.ngrid) - return gausslegendre_info(lobatto,radau,mass_matrix,sparse(S_matrix),K_matrix,L_matrix,mass_matrix_lu,Qmat) + + return gausslegendre_info(lobatto,radau,mass_matrix,sparse(S_matrix),sparse(K_matrix),sparse(L_matrix),mass_matrix_lu,Qmat) end -function setup_gausslegendre_pseudospectral_lobatto(coord;init_YY=true) +function setup_gausslegendre_pseudospectral_lobatto(coord; collision_operator_dim=true) x, w = gausslobatto(coord.ngrid) Dmat = allocate_float(coord.ngrid, coord.ngrid) gausslobattolegendre_differentiation_matrix!(Dmat,x,coord.ngrid) @@ -125,37 +132,38 @@ function setup_gausslegendre_pseudospectral_lobatto(coord;init_YY=true) GaussLegendre_weak_product_matrix!(M0,coord.ngrid,x,w,"M0") M1 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(M1,coord.ngrid,x,w,"M1") - M2 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(M2,coord.ngrid,x,w,"M2") - S0 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(S0,coord.ngrid,x,w,"S0") - S1 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(S1,coord.ngrid,x,w,"S1") K0 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(K0,coord.ngrid,x,w,"K0") K1 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(K1,coord.ngrid,x,w,"K1") - K2 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(K2,coord.ngrid,x,w,"K2") P0 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(P0,coord.ngrid,x,w,"P0") - P1 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(P1,coord.ngrid,x,w,"P1") - P2 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(P2,coord.ngrid,x,w,"P2") D0 = allocate_float(coord.ngrid) #@. D0 = Dmat[1,:] # values at lower extreme of element GaussLegendre_derivative_vector!(D0,-1.0,coord.ngrid,x,w) - Y00 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y01 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y10 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y11 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y20 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y21 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y30 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y31 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - if init_YY + if collision_operator_dim + M2 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(M2,coord.ngrid,x,w,"M2") + S0 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(S0,coord.ngrid,x,w,"S0") + S1 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(S1,coord.ngrid,x,w,"S1") + K2 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(K2,coord.ngrid,x,w,"K2") + P1 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(P1,coord.ngrid,x,w,"P1") + P2 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(P2,coord.ngrid,x,w,"P2") + + Y00 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y01 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y10 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y11 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y20 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y21 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y30 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y31 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(Y00,coord.ngrid,x,w,"Y00") GaussLegendre_weak_product_matrix!(Y01,coord.ngrid,x,w,"Y01") GaussLegendre_weak_product_matrix!(Y10,coord.ngrid,x,w,"Y10") @@ -164,12 +172,28 @@ function setup_gausslegendre_pseudospectral_lobatto(coord;init_YY=true) GaussLegendre_weak_product_matrix!(Y21,coord.ngrid,x,w,"Y21") GaussLegendre_weak_product_matrix!(Y30,coord.ngrid,x,w,"Y30") GaussLegendre_weak_product_matrix!(Y31,coord.ngrid,x,w,"Y31") + else + M2 = allocate_float(0, 0) + S0 = allocate_float(0, 0) + S1 = allocate_float(0, 0) + K2 = allocate_float(0, 0) + P1 = allocate_float(0, 0) + P2 = allocate_float(0, 0) + + Y00 = allocate_float(0, 0, 0) + Y01 = allocate_float(0, 0, 0) + Y10 = allocate_float(0, 0, 0) + Y11 = allocate_float(0, 0, 0) + Y20 = allocate_float(0, 0, 0) + Y21 = allocate_float(0, 0, 0) + Y30 = allocate_float(0, 0, 0) + Y31 = allocate_float(0, 0, 0) end return gausslegendre_base_info(Dmat,M0,M1,M2,S0,S1, K0,K1,K2,P0,P1,P2,D0,Y00,Y01,Y10,Y11,Y20,Y21,Y30,Y31) end -function setup_gausslegendre_pseudospectral_radau(coord;init_YY=true) +function setup_gausslegendre_pseudospectral_radau(coord; collision_operator_dim=true) # Gauss-Radau points on [-1,1) x, w = gaussradau(coord.ngrid) # Gauss-Radau points on (-1,1] @@ -182,35 +206,36 @@ function setup_gausslegendre_pseudospectral_radau(coord;init_YY=true) GaussLegendre_weak_product_matrix!(M0,coord.ngrid,xreverse,wreverse,"M0",radau=true) M1 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(M1,coord.ngrid,xreverse,wreverse,"M1",radau=true) - M2 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(M2,coord.ngrid,xreverse,wreverse,"M2",radau=true) - S0 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(S0,coord.ngrid,xreverse,wreverse,"S0",radau=true) - S1 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(S1,coord.ngrid,xreverse,wreverse,"S1",radau=true) K0 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(K0,coord.ngrid,xreverse,wreverse,"K0",radau=true) K1 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(K1,coord.ngrid,xreverse,wreverse,"K1",radau=true) - K2 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(K2,coord.ngrid,xreverse,wreverse,"K2",radau=true) P0 = allocate_float(coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(P0,coord.ngrid,xreverse,wreverse,"P0",radau=true) - P1 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(P1,coord.ngrid,xreverse,wreverse,"P1",radau=true) - P2 = allocate_float(coord.ngrid, coord.ngrid) - GaussLegendre_weak_product_matrix!(P2,coord.ngrid,xreverse,wreverse,"P2",radau=true) D0 = allocate_float(coord.ngrid) GaussLegendre_derivative_vector!(D0,-1.0,coord.ngrid,xreverse,wreverse,radau=true) - Y00 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y01 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y10 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y11 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y20 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y21 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y30 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - Y31 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) - if init_YY + if collision_operator_dim + M2 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(M2,coord.ngrid,xreverse,wreverse,"M2",radau=true) + S0 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(S0,coord.ngrid,xreverse,wreverse,"S0",radau=true) + S1 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(S1,coord.ngrid,xreverse,wreverse,"S1",radau=true) + K2 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(K2,coord.ngrid,xreverse,wreverse,"K2",radau=true) + P1 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(P1,coord.ngrid,xreverse,wreverse,"P1",radau=true) + P2 = allocate_float(coord.ngrid, coord.ngrid) + GaussLegendre_weak_product_matrix!(P2,coord.ngrid,xreverse,wreverse,"P2",radau=true) + + Y00 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y01 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y10 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y11 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y20 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y21 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y30 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) + Y31 = allocate_float(coord.ngrid, coord.ngrid, coord.ngrid) GaussLegendre_weak_product_matrix!(Y00,coord.ngrid,xreverse,wreverse,"Y00",radau=true) GaussLegendre_weak_product_matrix!(Y01,coord.ngrid,xreverse,wreverse,"Y01",radau=true) GaussLegendre_weak_product_matrix!(Y10,coord.ngrid,xreverse,wreverse,"Y10",radau=true) @@ -219,6 +244,22 @@ function setup_gausslegendre_pseudospectral_radau(coord;init_YY=true) GaussLegendre_weak_product_matrix!(Y21,coord.ngrid,xreverse,wreverse,"Y21",radau=true) GaussLegendre_weak_product_matrix!(Y30,coord.ngrid,xreverse,wreverse,"Y30",radau=true) GaussLegendre_weak_product_matrix!(Y31,coord.ngrid,xreverse,wreverse,"Y31",radau=true) + else + M2 = allocate_float(0, 0) + S0 = allocate_float(0, 0) + S1 = allocate_float(0, 0) + K2 = allocate_float(0, 0) + P1 = allocate_float(0, 0) + P2 = allocate_float(0, 0) + + Y00 = allocate_float(0, 0, 0) + Y01 = allocate_float(0, 0, 0) + Y10 = allocate_float(0, 0, 0) + Y11 = allocate_float(0, 0, 0) + Y20 = allocate_float(0, 0, 0) + Y21 = allocate_float(0, 0, 0) + Y30 = allocate_float(0, 0, 0) + Y31 = allocate_float(0, 0, 0) end return gausslegendre_base_info(Dmat,M0,M1,M2,S0,S1, K0,K1,K2,P0,P1,P2,D0,Y00,Y01,Y10,Y11,Y20,Y21,Y30,Y31) @@ -237,7 +278,7 @@ function elementwise_derivative!(coord, ff, gausslegendre::gausslegendre_info) imin = coord.imin[j]-k # imax is the maximum index on the full grid for this (jth) element imax = coord.imax[j] - if coord.name == "vperp" && coord.irank == 0 # differentiate this element with the Radau scheme + if coord.radau_first_element && coord.irank == 0 # differentiate this element with the Radau scheme @views mul!(df[:,j],gausslegendre.radau.Dmat[:,:],ff[imin:imax]) else #differentiate using the Lobatto scheme @views mul!(df[:,j],gausslegendre.lobatto.Dmat[:,:],ff[imin:imax]) @@ -267,81 +308,26 @@ function elementwise_derivative!(coord, ff, adv_fac, spectral::gausslegendre_inf return elementwise_derivative!(coord, ff, spectral) end -function elementwise_apply_Kmat!(coord, ff, gausslegendre::gausslegendre_info) - df = coord.scratch_2d - # define local variable nelement for convenience - nelement = coord.nelement_local - # check array bounds - @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df)) - - # variable k will be used to avoid double counting of overlapping point - k = 0 - j = 1 # the first element - imin = coord.imin[j]-k - # imax is the maximum index on the full grid for this (jth) element - imax = coord.imax[j] - get_KK_local!(gausslegendre.Qmat,j,gausslegendre.lobatto,gausslegendre.radau,coord,explicit_BC_terms=true) - #println(gausslegendre.Qmat) - @views mul!(df[:,j],gausslegendre.Qmat[:,:],ff[imin:imax]) - zero_gradient_bc_lower_boundary = false#true - if coord.name == "vperp" && zero_gradient_bc_lower_boundary - # set the 1st point of the RHS vector to zero - # consistent with use with the mass matrix with D f = 0 boundary conditions - df[1,j] = 0.0 - end - # calculate the derivative on each element - @inbounds for j ∈ 2:nelement - k = 1 - imin = coord.imin[j]-k - # imax is the maximum index on the full grid for this (jth) element - imax = coord.imax[j] - #@views mul!(df[:,j],gausslegendre.lobatto.Kmat[:,:],ff[imin:imax]) - get_KK_local!(gausslegendre.Qmat,j,gausslegendre.lobatto,gausslegendre.radau,coord,explicit_BC_terms=true) - #println(gausslegendre.Qmat) - @views mul!(df[:,j],gausslegendre.Qmat[:,:],ff[imin:imax]) - end - #for j in 1:nelement - # println(df[:,j]) - #end - return nothing -end +function single_element_interpolate!(result, newgrid, f, imin, imax, ielement, coord, + gausslegendre::gausslegendre_base_info) + n_new = length(newgrid) -function elementwise_apply_Lmat!(coord, ff, gausslegendre::gausslegendre_info) - df = coord.scratch_2d - # define local variable nelement for convenience - nelement = coord.nelement_local - # check array bounds - @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df)) - - # variable k will be used to avoid double counting of overlapping point - k = 0 - j = 1 # the first element - imin = coord.imin[j]-k - # imax is the maximum index on the full grid for this (jth) element - imax = coord.imax[j] - get_LL_local!(gausslegendre.Qmat,j,gausslegendre.lobatto,gausslegendre.radau,coord,explicit_BC_terms=true) - #println(gausslegendre.Qmat) - @views mul!(df[:,j],gausslegendre.Qmat[:,:],ff[imin:imax]) - zero_gradient_bc_lower_boundary = false#true - if coord.name == "vperp" && zero_gradient_bc_lower_boundary - # set the 1st point of the RHS vector to zero - # consistent with use with the mass matrix with D f = 0 boundary conditions - df[1,j] = 0.0 + i = 1 + other_nodes = @view coord.other_nodes[:,i,ielement] + one_over_denominator = coord.one_over_denominator[i,ielement] + this_f = f[i] + for j ∈ 1:n_new + result[j] = this_f * lagrange_poly_optimised(other_nodes, one_over_denominator, newgrid[j]) end - # calculate the derivative on each element - @inbounds for j ∈ 2:nelement - k = 1 - imin = coord.imin[j]-k - # imax is the maximum index on the full grid for this (jth) element - imax = coord.imax[j] - #@views mul!(df[:,j],gausslegendre.lobatto.Kmat[:,:],ff[imin:imax]) - get_LL_local!(gausslegendre.Qmat,j,gausslegendre.lobatto,gausslegendre.radau,coord,explicit_BC_terms=true) - #println(gausslegendre.Qmat) - @views mul!(df[:,j],gausslegendre.Qmat[:,:],ff[imin:imax]) + for i ∈ 2:coord.ngrid + other_nodes = @view coord.other_nodes[:,i,ielement] + one_over_denominator = coord.one_over_denominator[i,ielement] + this_f = f[i] + for j ∈ 1:n_new + result[j] += this_f * lagrange_poly_optimised(other_nodes, one_over_denominator, newgrid[j]) + end end - #for j in 1:nelement - # println(df[:,j]) - #end + return nothing end @@ -828,8 +814,7 @@ end """ A function that assigns the local weak-form matrices to a global array QQ_global for later solving weak form of required -1D equation. This function only supports fully local grids -that have coord.nelement_local = coord.nelement_global. +1D equation. The 'option' variable is a flag for choosing the type of matrix to be constructed. @@ -864,11 +849,16 @@ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2}, # N.B. QQ varies with ielement for vperp, but not vpa # a radau element is used for the vperp grid (see get_QQ_local!()) get_QQ_local!(QQ_j,j,lobatto,radau,coord,option) - QQ_global[imin[j],imin[j]:imax[j]] .+= QQ_j[1,:] + if coord.bc == "periodic" && coord.nrank == 1 + QQ_global[imax[end], imin[j]:imax[j]] .+= QQ_j[1,:] ./ 2.0 + QQ_global[imin[j],imin[j]:imax[j]] .+= QQ_j[1,:] ./ 2.0 + else + QQ_global[imin[j],imin[j]:imax[j]] .+= QQ_j[1,:] + end for k in 2:imax[j]-imin[j] QQ_global[k,imin[j]:imax[j]] .+= QQ_j[k,:] end - if coord.nelement_local > 1 + if coord.nelement_local > 1 || (coord.bc == "periodic" && coord.nrank == 1) QQ_global[imax[j],imin[j]:imax[j]] .+= QQ_j[ngrid,:]./2.0 else QQ_global[imax[j],imin[j]:imax[j]] .+= QQ_j[ngrid,:] @@ -883,7 +873,12 @@ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2}, end # upper boundary assembly on element if j == coord.nelement_local - QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:] + if coord.bc == "periodic" && coord.nrank == 1 + QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:] / 2.0 + QQ_global[imin[1],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:] / 2.0 + else + QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:] + end else QQ_global[imax[j],imin[j]-1:imax[j]] .+= QQ_j[ngrid,:]./2.0 end @@ -987,30 +982,32 @@ function get_KK_local!(QQ,ielement, if ielement > 1 || coord.irank > 0 # lobatto points @. QQ = (shift_factor/scale_factor)*lobatto.K0 + lobatto.K1 - lobatto.P0 # boundary terms from integration by parts - if explicit_BC_terms && ielement == 1 + if explicit_BC_terms && ielement == 1 && coord.irank == 0 imin = coord.imin[ielement] - 1 @. QQ[1,:] -= coord.grid[imin]*lobatto.Dmat[1,:]/scale_factor end - if explicit_BC_terms && ielement == nelement + if explicit_BC_terms && ielement == nelement && coord.irank == coord.nrank - 1 imax = coord.imax[ielement] @. QQ[coord.ngrid,:] += coord.grid[imax]*lobatto.Dmat[coord.ngrid,:]/scale_factor end else # radau points @. QQ = (shift_factor/scale_factor)*radau.K0 + radau.K1 - radau.P0 # boundary terms from integration by parts - if explicit_BC_terms && ielement == nelement + if explicit_BC_terms && ielement == nelement && coord.irank == coord.nrank - 1 imax = coord.imax[ielement] @. QQ[coord.ngrid,:] += coord.grid[imax]*radau.Dmat[coord.ngrid,:]/scale_factor end end else # assume integrals of form int^infty_-infty (.) d vpa @. QQ = lobatto.K0/scale_factor - # boundary terms from integration by parts - if explicit_BC_terms && ielement == 1 - @. QQ[1,:] -= lobatto.Dmat[1,:]/scale_factor - end - if explicit_BC_terms && ielement == nelement - @. QQ[coord.ngrid,:] += lobatto.Dmat[coord.ngrid,:]/scale_factor + if coord.bc !== "periodic" + # boundary terms from integration by parts + if explicit_BC_terms && ielement == 1 && coord.irank == 0 + @. QQ[1,:] -= lobatto.Dmat[1,:]/scale_factor + end + if explicit_BC_terms && ielement == nelement && coord.irank == coord.nrank - 1 + @. QQ[coord.ngrid,:] += lobatto.Dmat[coord.ngrid,:]/scale_factor + end end end return nothing @@ -1055,18 +1052,18 @@ function get_LL_local!(QQ,ielement, if ielement > 1 || coord.irank > 0 # lobatto points @. QQ = (shift_factor/scale_factor)*lobatto.K0 + lobatto.K1 # boundary terms from integration by parts - if explicit_BC_terms && ielement == 1 + if explicit_BC_terms && ielement == 1 && coord.irank == 0 imin = coord.imin[ielement] - 1 @. QQ[1,:] -= coord.grid[imin]*lobatto.Dmat[1,:]/scale_factor end - if explicit_BC_terms && ielement == nelement + if explicit_BC_terms && ielement == nelement && coord.irank == coord.nrank - 1 imax = coord.imax[ielement] @. QQ[coord.ngrid,:] += coord.grid[imax]*lobatto.Dmat[coord.ngrid,:]/scale_factor end else # radau points @. QQ = (shift_factor/scale_factor)*radau.K0 + radau.K1 # boundary terms from integration by parts - if explicit_BC_terms && ielement == nelement + if explicit_BC_terms && ielement == nelement && coord.irank == coord.nrank - 1 imax = coord.imax[ielement] @. QQ[coord.ngrid,:] += coord.grid[imax]*radau.Dmat[coord.ngrid,:]/scale_factor end @@ -1074,10 +1071,10 @@ function get_LL_local!(QQ,ielement, else # d^2 (.) d vpa^2 -- assume integrals of form int^infty_-infty (.) d vpa @. QQ = lobatto.K0/scale_factor # boundary terms from integration by parts - if explicit_BC_terms && ielement == 1 + if explicit_BC_terms && ielement == 1 && coord.irank == 0 @. QQ[1,:] -= lobatto.Dmat[1,:]/scale_factor end - if explicit_BC_terms && ielement == nelement + if explicit_BC_terms && ielement == nelement && coord.irank == coord.nrank - 1 @. QQ[coord.ngrid,:] += lobatto.Dmat[coord.ngrid,:]/scale_factor end end diff --git a/moment_kinetics/src/gyroaverages.jl b/moment_kinetics/src/gyroaverages.jl index f76edc742..b288e0da6 100644 --- a/moment_kinetics/src/gyroaverages.jl +++ b/moment_kinetics/src/gyroaverages.jl @@ -12,6 +12,7 @@ export gyroaverage_pdf! using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float, allocate_shared_float using ..array_allocation: allocate_int, allocate_shared_int +using ..lagrange_polynomials: lagrange_poly using ..looping using ..communication: MPISharedArray, comm_block, _block_synchronize @@ -246,30 +247,6 @@ function elementlist!(elist,coordlist,coord) return end -""" -Copy of function in fokker_planck_calculus.jl -Lagrange polynomial -args: -j - index of l_j from list of nodes -x_nodes - array of x node values -x - point where interpolated value is returned -""" -function lagrange_poly(j,x_nodes,x) - # get number of nodes - n = size(x_nodes,1) - # location where l(x0) = 1 - x0 = x_nodes[j] - # evaluate polynomial - poly = 1.0 - for i in 1:j-1 - poly *= (x - x_nodes[i])/(x0 - x_nodes[i]) - end - for i in j+1:n - poly *= (x - x_nodes[i])/(x0 - x_nodes[i]) - end - return poly -end - """ function for gyroaveraging a field of shape (z,r) and filling the result into an array of shape (vperp,z,r,s) diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 7f724c6bc..da0841fa3 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -636,6 +636,15 @@ function init_ion_pdf_over_density!(pdf, spec, composition, vpa, vperp, z, right_weight*upper_z_pdf_buffer) end + # Add a non-flowing Maxwellian (that vanishes at the sheath entrance boundaries) to try to + # avoid the 'hole' in the distribution function that can drive instabilities. + @loop_z_vperp iz ivperp begin + @. pdf[:,ivperp,iz] += spec.z_IC.density_amplitude * + (1.0 - (2.0 * z.grid[iz] / z.L)^2) * + exp(-(vpa.grid^2 + vperp.grid[ivperp]^2) + / vth[iz]^2) / vth[iz] + end + # Get the unnormalised pdf and the moments of the constructed full-f # distribution function (which will be modified from the input moments). convert_full_f_ion_to_normalised!(pdf, density, upar, ppar, vth, vperp, diff --git a/moment_kinetics/src/interpolation.jl b/moment_kinetics/src/interpolation.jl index 1e679ee10..af1ae1514 100644 --- a/moment_kinetics/src/interpolation.jl +++ b/moment_kinetics/src/interpolation.jl @@ -9,7 +9,24 @@ export interpolate_to_grid_z using ..array_allocation: allocate_float using ..moment_kinetics_structs: null_spatial_dimension_info, null_velocity_dimension_info -using ..type_definitions: mk_float +using ..type_definitions: mk_float, mk_int + +""" + single_element_interpolate!(result, newgrid, f, imin, imax, coord, spectral) + +Interpolation within a single element. + +`f` is an array with the values of the input variable in the element to be interpolated. +`imin` and `imax` give the start and end points of the element in the grid (used to +calculate shift and scale factors to a normalised grid). + +`newgrid` gives the points within the element where output is required. `result` is filled +with the interpolated values at those points. + +`coord` is the `coordinate` struct for the dimension along which interpolation is being +done. `spectral` is the corresponding `discretization_info`. +""" +function single_element_interpolate! end """ Interpolation from a regular grid to a 1d grid with arbitrary spacing @@ -30,6 +47,80 @@ spectral : discretization_info """ function interpolate_to_grid_1d! end +function interpolate_to_grid_1d!(result, newgrid, f, coord, spectral) + # define local variable nelement for convenience + nelement = coord.nelement_local + + n_new = size(newgrid)[1] + # Find which points belong to which element. + # kstart[j] contains the index of the first point in newgrid that is within element + # j, and kstart[nelement+1] is n_new if the last point is within coord.grid, or the + # index of the first element outside coord.grid otherwise. + # Assumes points in newgrid are sorted. + # May not be the most efficient algorithm. + # Find start/end points for each element, storing their indices in kstart + kstart = Vector{mk_int}(undef, nelement+1) + + # First element includes both boundary points, while all others have only one (to + # avoid duplication), so calculate the first element outside the loop. + if coord.radau_first_element && coord.irank == 0 + first_element_spectral = spectral.radau + # For a grid with a Radau first element, the lower boundary of the first element + # is at coord=0, and the coordinate range is 0 0.0 # centred second derivative for dissipation - @views derivative_z!(dummy_zrs, density, buffer_r_1, buffer_r_2, buffer_r_3, - buffer_r_4, z_spectral, z) - @views derivative_z!(moments.ion.d2dens_dz2, dummy_zrs, buffer_r_1, - buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) + @views second_derivative_z!(moments.ion.d2dens_dz2, density, buffer_r_1, + buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) end if moments.evolve_density || moments.evolve_upar || moments.evolve_ppar @views derivative_z!(moments.ion.dupar_dz, upar, buffer_r_1, @@ -853,12 +851,9 @@ function calculate_ion_moment_derivatives!(moments, scratch, scratch_dummy, z, z buffer_r_5, buffer_r_6, z_spectral, z) end if moments.evolve_upar && ion_mom_diss_coeff > 0.0 - # centred second derivative for dissipation - @views derivative_z!(dummy_zrs, upar, buffer_r_1, buffer_r_2, buffer_r_3, - buffer_r_4, z_spectral, z) - @views derivative_z!(moments.ion.d2upar_dz2, dummy_zrs, buffer_r_1, - buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) + @views second_derivative_z!(moments.ion.d2upar_dz2, upar, buffer_r_1, buffer_r_2, + buffer_r_3, buffer_r_4, z_spectral, z) end if moments.evolve_upar @views derivative_z!(moments.ion.dppar_dz, ppar, buffer_r_1, @@ -873,11 +868,11 @@ function calculate_ion_moment_derivatives!(moments, scratch, scratch_dummy, z, z buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4, buffer_r_5, buffer_r_6, z_spectral, z) - # centred second derivative for dissipation - @views derivative_z!(dummy_zrs, ppar, buffer_r_1, buffer_r_2, buffer_r_3, - buffer_r_4, z_spectral, z) - @views derivative_z!(moments.ion.d2ppar_dz2, dummy_zrs, buffer_r_1, - buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) + if ion_mom_diss_coeff > 0.0 + # centred second derivative for dissipation + @views second_derivative_z!(moments.ion.d2ppar_dz2, ppar, buffer_r_1, + buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) + end @views derivative_z!(moments.ion.dqpar_dz, qpar, buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) @@ -1391,11 +1386,9 @@ function calculate_neutral_moment_derivatives!(moments, scratch, scratch_dummy, if moments.evolve_density && neutral_mom_diss_coeff > 0.0 # centred second derivative for dissipation - @views derivative_z!(dummy_zrsn, density, buffer_r_1, buffer_r_2, buffer_r_3, - buffer_r_4, z_spectral, z; neutrals=true) - @views derivative_z!(moments.neutral.d2dens_dz2, dummy_zrsn, - buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4, - z_spectral, z; neutrals=true) + @views second_derivative_z!(moments.neutral.d2dens_dz2, density, buffer_r_1, + buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z; + neutrals=true) end if moments.evolve_density || moments.evolve_upar || moments.evolve_ppar @views derivative_z!(moments.neutral.duz_dz, uz, buffer_r_1, @@ -1413,12 +1406,9 @@ function calculate_neutral_moment_derivatives!(moments, scratch, scratch_dummy, buffer_r_5, buffer_r_6, z_spectral, z; neutrals=true) end if moments.evolve_upar && neutral_mom_diss_coeff > 0.0 - # centred second derivative for dissipation - @views derivative_z!(dummy_zrsn, uz, buffer_r_1, buffer_r_2, buffer_r_3, - buffer_r_4, z_spectral, z; neutrals=true) - @views derivative_z!(moments.neutral.d2uz_dz2, dummy_zrsn, buffer_r_1, - buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z) + @views second_derivative_z!(moments.neutral.d2uz_dz2, uz, buffer_r_1, buffer_r_2, + buffer_r_3, buffer_r_4, z_spectral, z; neutrals=true) end if moments.evolve_upar @views derivative_z!(moments.neutral.dpz_dz, pz, buffer_r_1, @@ -1434,12 +1424,12 @@ function calculate_neutral_moment_derivatives!(moments, scratch, scratch_dummy, buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4, buffer_r_5, buffer_r_6, z_spectral, z; neutrals=true) - # centred second derivative for dissipation - @views derivative_z!(dummy_zrsn, pz, buffer_r_1, buffer_r_2, buffer_r_3, - buffer_r_4, z_spectral, z; neutrals=true) - @views derivative_z!(moments.neutral.d2pz_dz2, dummy_zrsn, buffer_r_1, - buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z; - neutrals=true) + if neutral_mom_diss_coeff > 0.0 + # centred second derivative for dissipation + @views second_derivative_z!(moments.neutral.d2pz_dz2, pz, buffer_r_1, + buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z; + neutrals=true) + end @views derivative_z!(moments.neutral.dqz_dz, qz, buffer_r_1, buffer_r_2, buffer_r_3, buffer_r_4, z_spectral, z; diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl index 52c0aee90..e439ac4cf 100644 --- a/moment_kinetics/test/calculus_tests.jl +++ b/moment_kinetics/test/calculus_tests.jl @@ -904,7 +904,7 @@ function runtests() # create the coordinate struct 'x' and info for derivatives, etc. # This test runs effectively in serial, so use `ignore_MPI=true` to avoid # errors due to communicators not being fully set up. - x, spectral = define_coordinate(input; ignore_MPI=true, init_YY=false) + x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -1023,7 +1023,7 @@ function runtests() # create the coordinate struct 'x' and info for derivatives, etc. # This test runs effectively in serial, so use `ignore_MPI=true` to avoid # errors due to communicators not being fully set up. - x, spectral = define_coordinate(input; ignore_MPI=true, init_YY=false) + x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset @@ -1069,7 +1069,7 @@ function runtests() # create the coordinate struct 'x' and info for derivatives, etc. # This test runs effectively in serial, so use `ignore_MPI=true` to avoid # errors due to communicators not being fully set up. - x, spectral = define_coordinate(input; ignore_MPI=true, init_YY=false) + x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) # test polynomials up to order ngrid-1 for n ∈ 0:ngrid-1 # create array for the function f(x) to be differentiated/integrated @@ -1122,7 +1122,7 @@ function runtests() # create the coordinate struct 'x' and info for derivatives, etc. # This test runs effectively in serial, so use `ignore_MPI=true` to avoid # errors due to communicators not being fully set up. - x, spectral = define_coordinate(input; ignore_MPI=true, init_YY=false) + x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) # test polynomials up to order ngrid-1 for n ∈ 0:ngrid-1 # create array for the function f(x) to be differentiated/integrated @@ -1448,7 +1448,7 @@ function runtests() # create the coordinate struct 'x' and info for derivatives, etc. # This test runs effectively in serial, so use `ignore_MPI=true` to avoid # errors due to communicators not being fully set up. - x, spectral = define_coordinate(input; ignore_MPI=true, init_YY=false) + x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) offset = randn(rng) f = @. sinpi(2.0 * x.grid / L) + offset diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl index 3039b934d..08239c0f7 100644 --- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl +++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl @@ -39,47 +39,47 @@ const expected = [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], # Expected phi: - [-1.267505494648937, -1.275683298550937], + [-1.268504994982021, -1.276675261324553], # Expected n_ion: - [0.2815330322340072, 0.2792400986636072], + [0.28125178045355353, 0.278963240220169], # Expected upar_ion: [0.0, 0.0], # Expected ppar_ion: - [0.17982280248048935, 0.14891126175332367], + [0.17964315932116812, 0.148762861611732], # Expected pperp_ion - [0.14340146667506784, 0.1581377822859991], + [0.14325820846660123, 0.15798027481288696], # Expected qpar_ion [0.0, 0.0], # Expected v_t_ion - [1.0511726083010418, 1.0538509291794658], + [1.0511726083010418, 1.0538484394097123], # Expected dSdt - [0.0, 1.1853081348031516e-5], + [0.0, 1.1831920390587679e-5], # Expected f_ion: [0.0 0.0 0.0 0.0 0.0 0.0 0.0; - 0.0006199600161806666 0.00047805300997075977 0.0002665817112117718 7.637693901737056e-5 1.3272321881722645e-5 1.3988924344690309e-6 0.0; - 0.005882016862626724 0.0045356406743786385 0.002529256854781707 0.0007246442213864763 0.00012592428394890537 1.3272321881722645e-5 0.0; - 0.03384866997225574 0.026100809957763767 0.01455486826237011 0.004170039574837177 0.000724644221263874 7.637693900444835e-5 0.0; - 0.11813810317200342 0.09109664226661075 0.05079917556123135 0.01455420747483572 0.0025291420266981487 0.0002665696084208068 0.0; - 0.22957946520936198 0.17702940755267918 0.0987187642706944 0.0282833995036643 0.004914917865936108 0.0005180285318549189 0.0; + 0.0006193406755051616 0.0004775754345362236 0.0002663153958159559 7.630063837899157e-5 1.3259062818903743e-5 1.3974949395295016e-6 0.0; + 0.005876140721904817 0.0045311095648138235 0.00252673012465705 0.000723920301085391 0.00012579848546344195 1.3259062818903743e-5 0.0; + 0.0338148551171386 0.026074735222541223 0.01454032793443568 0.004165873701136042 0.000723920300962911 7.630063836608227e-5 0.0; + 0.11802008308891455 0.09100563662998079 0.05074842713409726 0.014539667807028695 0.0025266154112868616 0.0002663033051156912 0.0; + 0.22935011509426767 0.1768525549976815 0.09862014412656786 0.028255144359305 0.00491000785807803 0.0005175110208340848 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; - 0.22957946520936204 0.1770294075526792 0.0987187642706944 0.0282833995036643 0.004914917865936108 0.0005180285318549189 0.0; - 0.11813810317200349 0.0910966422666108 0.050799175561231376 0.01455420747483573 0.0025291420266981487 0.0002665696084208068 0.0; - 0.03384866997225574 0.026100809957763767 0.01455486826237011 0.004170039574837177 0.000724644221263874 7.637693900444835e-5 0.0; - 0.005882016862626724 0.0045356406743786385 0.002529256854781707 0.0007246442213864763 0.00012592428394890537 1.3272321881722645e-5 0.0; - 0.0006199600161806666 0.00047805300997075977 0.0002665817112117718 7.637693901737056e-5 1.3272321881722645e-5 1.3988924344690309e-6 0.0; + 0.22935011509426778 0.17685255499768154 0.09862014412656786 0.028255144359305 0.00491000785807803 0.0005175110208340848 0.0; + 0.1180200830889146 0.09100563662998083 0.05074842713409728 0.014539667807028703 0.0025266154112868616 0.0002663033051156912 0.0; + 0.0338148551171386 0.026074735222541223 0.01454032793443568 0.004165873701136042 0.000723920300962911 7.630063836608227e-5 0.0; + 0.005876140721904817 0.0045311095648138235 0.00252673012465705 0.000723920301085391 0.00012579848546344195 1.3259062818903743e-5 0.0; + 0.0006193406755051616 0.0004775754345362236 0.0002663153958159559 7.630063837899157e-5 1.3259062818903743e-5 1.3974949395295016e-6 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0;;; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; - 0.0001712743622973216 7.105465094508053e-5 -7.829380680167827e-5 -0.00015364081956318698 -9.097098213067502e-5 -3.311284120491419e-5 0.0; - 0.005883280697248667 0.004667594200766182 0.002855965521103658 0.0008138347136178689 2.44260649525292e-5 -9.753249634264602e-5 0.0; - 0.02792209301450194 0.022385716644538384 0.01413535091105969 0.004677801530322722 0.0007105315221401102 -0.00022400635166536323 0.0; - 0.08117458037332098 0.06563459159004267 0.04247673844050208 0.015087784332275832 0.0029056314178876035 -0.00023019804543218203 0.0; - 0.15133793170654106 0.12313903060106579 0.08111673445361306 0.029975277983613262 0.00626735398468981 7.553501812465833e-6 0.0; - 0.18493902160817713 0.15073513412904313 0.09976414473955808 0.037251581926306565 0.007941836186495122 0.00016196175024033304 0.0; - 0.15133793170654092 0.12313903060106571 0.08111673445361306 0.02997527798361324 0.006267353984689816 7.553501812469816e-6 0.0; - 0.081174580373321 0.06563459159004267 0.042476738440502065 0.015087784332275821 0.002905631417887614 -0.0002301980454321778 0.0; - 0.027922093014501933 0.022385716644538384 0.014135350911059698 0.004677801530322729 0.0007105315221401184 -0.00022400635166536134 0.0; - 0.005883280697248667 0.004667594200766184 0.002855965521103663 0.0008138347136178759 2.4426064952530956e-5 -9.753249634264635e-5 0.0; - 0.0001712743622973275 7.105465094508572e-5 -7.829380680167411e-5 -0.00015364081956318568 -9.097098213067551e-5 -3.311284120491447e-5 0.0; + 0.00017108987342944037 7.097261590252227e-5 -7.822316408658004e-5 -0.000153489754637506 -9.087984332447761e-5 -3.307937957312587e-5 0.0; + 0.005877384425022921 0.004662912823128315 0.002853094592035005 0.0008130106752860298 2.4399432485772724e-5 -9.74348090421241e-5 0.0; + 0.02789429969714529 0.022363413859309983 0.014121230173998248 0.004673094872084004 0.0007098078939540657 -0.0002237857510708618 0.0; + 0.08109427896718696 0.06556961121847364 0.04243458963422585 0.01507272286376149 0.0029027058292680945 -0.0002299783231637593 0.0; + 0.1511887162169901 0.12301753182687214 0.08103653941734254 0.029945484317773212 0.0062610744742386155 7.528837370841561e-6 0.0; + 0.184756848099325 0.1505865499710698 0.09966561578213134 0.037214599991686845 0.007933889035324836 0.00016178010211991204 0.0; + 0.1511887162169905 0.12301753182687247 0.08103653941734275 0.029945484317773295 0.006261074474238628 7.528837370846196e-6 0.0; + 0.08109427896718732 0.06556961121847393 0.042434589634226055 0.015072722863761552 0.0029027058292681127 -0.0002299783231637549 0.0; + 0.027894299697145436 0.022363413859310108 0.014121230173998337 0.00467309487208404 0.0007098078939540783 -0.00022378575107086105 0.0; + 0.005877384425022965 0.00466291282312835 0.0028530945920350282 0.0008130106752860425 2.4399432485774198e-5 -9.743480904212476e-5 0.0; + 0.00017108987342944414 7.097261590252536e-5 -7.822316408657793e-5 -0.0001534897546375058 -9.087984332447822e-5 -3.3079379573126077e-5 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0]) ########################################################################################### # to modify the test, with a new expected f, print the new f using the following commands @@ -135,7 +135,7 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral", "electron_physics" => "boltzmann_electron_response", "fokker_planck_collisions" => Dict{String,Any}("use_fokker_planck" => true, "nuii" => 1.0, "frequency_option" => "manual"), "z_IC_upar_amplitude1" => 0.0, - "z_IC_density_amplitude1" => 0.001, + "z_IC_density_amplitude1" => 0.0, "z_IC_upar_amplitude2" => 0.0, "z_IC_temperature_phase1" => 0.0, "z_IC_temperature_amplitude1" => 0.0, diff --git a/moment_kinetics/test/gyroaverage_tests.jl b/moment_kinetics/test/gyroaverage_tests.jl index 8b1740d02..38a6c7085 100644 --- a/moment_kinetics/test/gyroaverage_tests.jl +++ b/moment_kinetics/test/gyroaverage_tests.jl @@ -135,11 +135,11 @@ function gyroaverage_test(absolute_error; rhostar=0.1, pitch=0.5, ngrid=5, kr=2, nrank, irank, gyrophase_L, gyrophase_discretization, fd_option, cheb_option, "periodic", adv_input,comm,element_spacing_option) # create the coordinate structs - r, r_spectral = define_coordinate(r_input; init_YY=false, run_directory=test_output_directory) - z, z_spectral = define_coordinate(z_input; init_YY=false, run_directory=test_output_directory) - vperp, vperp_spectral = define_coordinate(vperp_input; init_YY=false, run_directory=test_output_directory) - vpa, vpa_spectral = define_coordinate(vpa_input; init_YY=false, run_directory=test_output_directory) - gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input; init_YY=false, run_directory=test_output_directory) + r, r_spectral = define_coordinate(r_input; collision_operator_dim=false, run_directory=test_output_directory) + z, z_spectral = define_coordinate(z_input; collision_operator_dim=false, run_directory=test_output_directory) + vperp, vperp_spectral = define_coordinate(vperp_input; collision_operator_dim=false, run_directory=test_output_directory) + vpa, vpa_spectral = define_coordinate(vpa_input; collision_operator_dim=false, run_directory=test_output_directory) + gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input; collision_operator_dim=false, run_directory=test_output_directory) # create test geometry #rhostar = 0.1 #rhostar of ions for ExB drift diff --git a/moment_kinetics/test/harrisonthompson.jl b/moment_kinetics/test/harrisonthompson.jl index 553515a19..4e4cada7d 100644 --- a/moment_kinetics/test/harrisonthompson.jl +++ b/moment_kinetics/test/harrisonthompson.jl @@ -74,7 +74,7 @@ test_input_finite_difference = Dict("n_ion_species" => 1, "initial_density1" => 1.0, "initial_temperature1" => 1.0, "z_IC_option1" => "gaussian", - "z_IC_density_amplitude1" => 0.001, + "z_IC_density_amplitude1" => 0.0, "z_IC_density_phase1" => 0.0, "z_IC_upar_amplitude1" => 0.0, "z_IC_upar_phase1" => 0.0, diff --git a/moment_kinetics/test/interpolation_tests.jl b/moment_kinetics/test/interpolation_tests.jl index 9b88deda7..5e9a1d59c 100644 --- a/moment_kinetics/test/interpolation_tests.jl +++ b/moment_kinetics/test/interpolation_tests.jl @@ -26,10 +26,12 @@ adv_input = advection_input("default", 1.0, 0.0, 0.0) function runtests() @testset "interpolation" verbose=use_verbose begin - @testset "$discretization, $ntest, $nelement, $zlim" for + @testset "$discretization, $element_spacing_option, $ntest, $nelement, $zlim" for (discretization, element_spacing_option, rtol) ∈ (("finite_difference", "uniform", 1.e-5), ("chebyshev_pseudospectral", "uniform", 1.e-8), - ("chebyshev_pseudospectral", "sqrt", 1.e-8)), + ("chebyshev_pseudospectral", "sqrt", 1.e-8), + ("gausslegendre_pseudospectral", "uniform", 1.e-8), + ("gausslegendre_pseudospectral", "sqrt", 1.e-8)), ntest ∈ (3, 14), nelement ∈ (2, 8), zlim ∈ (L/2.0, L/5.0) # create the 'input' struct containing input info needed to create a coordinate @@ -46,7 +48,7 @@ function runtests() # create the coordinate struct 'z' # This test runs effectively in serial, so use `ignore_MPI=true` to avoid # errors due to communicators not being fully set up. - z, spectral = define_coordinate(input; ignore_MPI=true) + z, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) test_grid = [z for z in range(-zlim, zlim, length=ntest)] diff --git a/moment_kinetics/test/recycling_fraction_tests.jl b/moment_kinetics/test/recycling_fraction_tests.jl index 9e3c1697d..935ed968e 100644 --- a/moment_kinetics/test/recycling_fraction_tests.jl +++ b/moment_kinetics/test/recycling_fraction_tests.jl @@ -32,7 +32,7 @@ test_input = Dict("n_ion_species" => 1, "initial_density1" => 1.0, "initial_temperature1" => 1.0, "z_IC_option1" => "gaussian", - "z_IC_density_amplitude1" => 0.001, + "z_IC_density_amplitude1" => 0.0, "z_IC_density_phase1" => 0.0, "z_IC_upar_amplitude1" => 1.0, "z_IC_upar_phase1" => 0.0, diff --git a/moment_kinetics/test/wall_bc_tests.jl b/moment_kinetics/test/wall_bc_tests.jl index 6c46e294e..c1c40634d 100644 --- a/moment_kinetics/test/wall_bc_tests.jl +++ b/moment_kinetics/test/wall_bc_tests.jl @@ -30,7 +30,7 @@ test_input_finite_difference = Dict("n_ion_species" => 1, "initial_density1" => 1.0, "initial_temperature1" => 1.0, "z_IC_option1" => "gaussian", - "z_IC_density_amplitude1" => 0.001, + "z_IC_density_amplitude1" => 0.0, "z_IC_density_phase1" => 0.0, "z_IC_upar_amplitude1" => 0.0, "z_IC_upar_phase1" => 0.0,