# Solving the Boltzmann transport equation using Trixi.jl

**URL:** <https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509>\
**Category:** Modelling & Simulations\
**Tags:** pde, trixi\
**Created:** [September 19, 2025, 10:06am UTC](https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509 "2025-09-19T10:06:45Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![drandran12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drandran12/32/220647_2.png) [@drandran12](https://discourse.julialang.org/u/drandran12)\
**Post date:** [September 19, 2025, 10:06am UTC](https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509/1 "2025-09-19T10:06:45Z")

</div>

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](https://github.com/ranocha/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?

```julia
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);

```

---

<div class="post-metadata">

**Author:** ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)\
**Post date:** [September 19, 2025, 11:31am UTC](https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509/2 "2025-09-19T11:31:33Z")

</div>

Trixi uses some kind of interface for the equations you can look at TrixiShallowWater.jl or BloodFlowTrixi.jl github code to see how specific models can be defined its all pure julia though so it should be compatible with everything you can think of, there is also tons of exemple in [Trixi.jl/src/equations at main · trixi-framework/Trixi.jl · GitHub](https://github.com/trixi-framework/Trixi.jl/tree/main/src/equations), the only question is will you have access to the variable/derivative you want at the mesh point you want which is not sure

---

<div class="post-metadata">

**Author:** ![JoshuaLampert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joshualampert/32/205964_2.png) [@JoshuaLampert](https://discourse.julialang.org/u/JoshuaLampert)\
**Post date:** [September 19, 2025, 12:11pm UTC](https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509/3 "2025-09-19T12:11:15Z")

</div>

A tutorial on defining your own equations, is given in [11 Adding a new scalar conservation law · Trixi.jl](https://trixi-framework.github.io/TrixiDocumentation/stable/tutorials/adding_new_scalar_equations/). The main part consists in defining the `flux`, which would be similar to your `advection_term!`. Defining this already allows you to use `flux_central` as numerical flux; if you want to use more sophisticated fluxes you would need to implement more. Your collision integral could then be implemented as a source term. The integral in there is a little bit tricky though.

---

<div class="post-metadata">

**Author:** ![drandran12](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/drandran12/32/220647_2.png) [@drandran12](https://discourse.julialang.org/u/drandran12)\
**Post date:** [September 19, 2025, 12:21pm UTC](https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509/4 "2025-09-19T12:21:35Z")

</div>

@JoshuaLampert Totally missed that tutorial. Greatly appreciate it.

I guess my main problem is that I have a “hidden” 4d problem: 2 spatial variables and associated velocity variables. Is there an interface for this? Should I somehow cast it into the `LatticeBoltzmannEquations2D` type equation?

---

<div class="post-metadata">

**Author:** ![JoshuaLampert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joshualampert/32/205964_2.png) [@JoshuaLampert](https://discourse.julialang.org/u/JoshuaLampert)\
**Post date:** [September 19, 2025, 12:35pm UTC](https://discourse.julialang.org/t/solving-the-boltzmann-transport-equation-using-trixi-jl/132509/5 "2025-09-19T12:35:58Z")

</div>

The 4d nature of the problem is an issue. True 4d is currently not possible with Trixi.jl. If Lattice-Boltzmann would be enough, you could try them instead. I don’t have much experience with these equations though.