I need to solve the two spatial dimensional Boltzmann transport equation of the form:
\left[\partial_t+\vec{k} \cdot \vec{\nabla}_x-\vec{\nabla}_x V \cdot \vec{\nabla}_k\right] n(\vec{x}, t, \vec{k})=-4 \Gamma_a \int \frac{d^d q}{(2 \pi)^d} n(\vec{x}, t, \vec{k}) n(\vec{x}, t, \vec{q})
with a specific collision integral on the right-hand-side. Using SummationByPartsOperators.jl I have implemented a naive implementation for the homogeneous case (see below). However, I would like to also attempt to make a implementation using Trixi.jl. Trixi tutorials seem to always use their prebuild equations. I find it hard to understand how to implement your own equation as you can do with SummationByPartsOperators.jl.
Can somebody point me in the right directions to get me started?
using SummationByPartsOperators, OrdinaryDiffEq, ProgressLogging
using LinearAlgebra: norm
using Plots: Plots, scatter, @animate, gif
using Printf
# Domain parameters
xmin, xmax = -2.0, 2.0
ymin, ymax = -2.0, 2.0
N_x, N_y = 32, 32
kxmin, kxmax = -π, π
kymin, kymax = -π, π
N_kx, N_ky = 16, 16
# Only need spatial derivatives - no k derivatives!
D_x = derivative_operator(
MattssonNordström2004();
derivative_order=1,
accuracy_order=4,
xmin=xmin,
xmax=xmax,
N=N_x,
)
D_y = derivative_operator(
MattssonNordström2004();
derivative_order=1,
accuracy_order=4,
xmin=ymin,
xmax=ymax,
N=N_y,
)
D_kx = derivative_operator(
MattssonNordström2004();
derivative_order=1,
accuracy_order=4,
xmin=kxmin,
xmax=kxmax,
N=N_kx,
)
D_ky = derivative_operator(
MattssonNordström2004();
derivative_order=1,
accuracy_order=4,
xmin=kymin,
xmax=kymax,
N=N_ky,
)
# Create spatial tensor product
D_spatial = tensor_product_operator_2D(D_x, D_y)
D_momentum = tensor_product_operator_2D(D_kx, D_ky)
# Grid points
nodes_spatial = SummationByPartsOperators.grid(D_spatial)
nodes_momentum = SummationByPartsOperators.grid(D_momentum)
N_total = N_x * N_y * N_kx * N_ky
x0(x) = exp(-10 * norm(x .- (0.0, 0.0))^2)
k0(k) = 1.0#exp(-30 * norm(k .- (1.0, 1.0))^2)
function initial_condition(x, k)
return x0(x) #* k0(k)
end
# Initialize solution
u0 = zeros(N_x * N_y, N_kx * N_ky)
for i in eachindex(nodes_spatial)
for j in eachindex(nodes_momentum)
u0[i, j] = initial_condition(nodes_spatial[i], nodes_momentum[j])
end
end
scatter(first.(nodes_spatial), last.(nodes_spatial), u0[:, 2]; zrange=(-0.1, 1.1))
# scatter(first.(nodes_momentum), last.(nodes_momentum), u0[1822, :]; zrange=(-0.1, 1.1))
# Collision integral - only integration, no derivatives
function collision_integral!(col, u, Γ_a, D_momentum)
for i in axes(u, 1) # iterate over space
density = integrate(u[i, :], D_momentum)
for j in axes(u, 2)
col[i, j] = -4 * Γ_a * density * u[i, j]
end
end
end
# Transport term - only spatial derivatives
function advection_term!(adv, u, D_spatial, nodes_momentum)
for i in axes(u, 2) # iterate over momentum
kv = nodes_momentum[i]
D = kv[1] * D_spatial[1] + kv[2] * D_spatial[2]
adv[:, i] .= -1.0 * D * u[:, i]
# mul!(du[:, i], D, u[:, i], -1.0, 0.0)
end
end
# Main ODE function
function boltzmann_2d!(du, u, p, t)
col, adv, Γ_a, D_momentum, D_spatial, nodes_momentum = p
collision_integral!(col, u, Γ_a, D_momentum)
advection_term!(adv, u, D_spatial, nodes_momentum)
@. du = col + adv
return nothing
end
# Solve
Γ_a = 0.1
col = zeros(N_x * N_y, N_kx * N_ky)
adv = zeros(N_x * N_y, N_kx * N_ky)
p = (col, adv, Γ_a, D_momentum, D_spatial, nodes_momentum)
tspan = (0.0, 0.5) # Shorter time for testing
prob = ODEProblem(boltzmann_2d!, u0, tspan, p)
frames = 50
saveat = range(tspan...; length=frames)
alg = Tsit5() # SSPRK104()
sol = solve(prob, alg; saveat, progress=true, progress_steps=100)
scatter(first.(nodes_spatial), last.(nodes_spatial), sol(0.0)[:, 120]; zrange=(-0.1, 1.1))
_, lowest_momentum_idx = findmin(norm, nodes_momentum)
anim = @animate for t in saveat
scatter(
first.(nodes_spatial),
last.(nodes_spatial),
sol(t)[:, lowest_momentum_idx];
title=@sprintf("t = %.2f", t),
zrange=(-0.1, 1.1),
label=@sprintf("k = (%.2f,%.2f)", nodes_momentum[lowest_momentum_idx]...),
)
end
gif(anim, "Boltzmann_2D.gif"; fps=frames Ă· 3);