Skip to content

Commit

Permalink
Change vperp numerical dissipation to use laplacian_derivative!(), su…
Browse files Browse the repository at this point in the history
…pport laplacian_derivative!() for chebyshev points with direct differentiation calculation.
mrhardman committed Jul 5, 2024
1 parent e6bdb90 commit 8c84a1e
Showing 2 changed files with 70 additions and 4 deletions.
70 changes: 68 additions & 2 deletions moment_kinetics/src/calculus.jl
Original file line number Diff line number Diff line change
@@ -131,6 +131,70 @@ function second_derivative!(d2f, f, coord, spectral; handle_periodic=true)
return nothing
end

function laplacian_derivative!(d2f, f, coord, spectral; handle_periodic=true)
# computes (1/coord) d / coord ( coord d f / d(coord)) for vperp coordinate
# For spectral element methods, calculate second derivative by applying first
# derivative twice, with special treatment for element boundaries

# First derivative
elementwise_derivative!(coord, f, spectral)
derivative_elements_to_full_grid!(coord.scratch3, coord.scratch_2d, coord)
if coord.name == "vperp"
# include the Jacobian
@. coord.scratch3 *= coord.grid
end
# MPI reconcile code here if used with z or r coords

# Save elementwise first derivative result
coord.scratch2_2d .= coord.scratch_2d

# Second derivative for element interiors
elementwise_derivative!(coord, coord.scratch3, spectral)
derivative_elements_to_full_grid!(d2f, coord.scratch_2d, coord)
if coord.name == "vperp"
# include the Jacobian
@. d2f /= coord.grid
end
# MPI reconcile code here if used with z or r coords

# Add contribution to penalise discontinuous first derivatives at element
# boundaries. For smooth functions this would do nothing so should not affect
# convergence of the second derivative. Aims to stabilise numerical instability when
# spike develops at an element boundary. The coefficient is an arbitrary choice, it
# should probably be large enough for stability but as small as possible.
#
# Arbitrary numerical coefficient
C = 1.0
function penalise_discontinuous_first_derivative!(d2f, imin, imax, df)
# Left element boundary
d2f[imin] += C * df[1]

# Right element boundary
d2f[imax] -= C * df[end]

return nothing
end
@views penalise_discontinuous_first_derivative!(d2f, 1, coord.imax[1],
coord.scratch2_2d[:,1])
for ielement 2:coord.nelement_local
@views penalise_discontinuous_first_derivative!(d2f, coord.imin[ielement]-1,
coord.imax[ielement],
coord.scratch2_2d[:,ielement])
end

if coord.bc ("zero", "both_zero", "zero-no-regularity")
# For stability don't contribute to evolution at boundaries, in case these
# points are not set by a boundary condition.
# Full grid may be across processes and bc only applied to extreme ends of the
# domain.
if coord.irank == coord.nrank - 1
d2f[end] = 0.0
end
else
error("Unsupported bc '$(coord.bc)'")
end
return nothing
end
"""
mass_matrix_solve!(f, b, spectral::weak_discretization_info)
@@ -176,9 +240,11 @@ function laplacian_derivative!(d2f, f, coord, spectral::weak_discretization_info
# for all other coord.name, do exactly the same as second_derivative! above.
mul!(coord.scratch, spectral.L_matrix, f)

if handle_periodic && coord.bc == "periodic"
if coord.bc == "periodic" && coord.name == "vperp"
error("laplacian_derivative!() cannot handle periodic boundaries for vperp")
elseif coord.bc == "periodic"
if coord.nrank > 1
error("second_derivative!() cannot handle periodic boundaries for a "
error("laplacian_derivative!() cannot handle periodic boundaries for a "
* "distributed coordinate")
end

4 changes: 2 additions & 2 deletions moment_kinetics/src/numerical_dissipation.jl
Original file line number Diff line number Diff line change
@@ -12,7 +12,7 @@ export setup_numerical_dissipation, vpa_boundary_buffer_decay!,
using Base.Iterators: flatten

using ..looping
using ..calculus: derivative!, second_derivative!
using ..calculus: derivative!, second_derivative!, laplacian_derivative!
using ..derivatives: derivative_r!, derivative_z!, second_derivative_r!,
second_derivative_z!
using ..type_definitions: mk_float
@@ -369,7 +369,7 @@ function vperp_dissipation!(f_out, f_in, vperp, spectral::T_spectral, dt,
end

@loop_s_r_z_vpa is ir iz ivpa begin
@views second_derivative!(vperp.scratch, f_in[ivpa,:,iz,ir,is], vperp, spectral)
@views laplacian_derivative!(vperp.scratch, f_in[ivpa,:,iz,ir,is], vperp, spectral)
@views @. f_out[ivpa,:,iz,ir,is] += dt * diffusion_coefficient * vperp.scratch
end

0 comments on commit 8c84a1e

Please sign in to comment.