diff --git a/src/equations/hyperbolic_diffusion_2d.jl b/src/equations/hyperbolic_diffusion_2d.jl
index 28f0c2bf47..a243f5b2a3 100644
--- a/src/equations/hyperbolic_diffusion_2d.jl
+++ b/src/equations/hyperbolic_diffusion_2d.jl
@@ -174,6 +174,12 @@ end
     sqrt(equations.nu * equations.inv_Tr) * norm(normal_direction)
 end
 
+"""
+    flux_godunov(u_ll, u_rr, orientation_or_normal_direction, 
+                 equations::HyperbolicDiffusionEquations2D)
+
+Godunov (upwind) flux for the 2D hyperbolic diffusion equations.
+"""
 @inline function flux_godunov(u_ll, u_rr, orientation::Integer,
                               equations::HyperbolicDiffusionEquations2D)
     # Obtain left and right fluxes
diff --git a/src/equations/hyperbolic_diffusion_3d.jl b/src/equations/hyperbolic_diffusion_3d.jl
index e9017b1544..030a0e1813 100644
--- a/src/equations/hyperbolic_diffusion_3d.jl
+++ b/src/equations/hyperbolic_diffusion_3d.jl
@@ -188,6 +188,12 @@ end
     λ_max = sqrt(equations.nu * equations.inv_Tr)
 end
 
+"""
+    flux_godunov(u_ll, u_rr, orientation_or_normal_direction, 
+                 equations::HyperbolicDiffusionEquations3D)
+
+Godunov (upwind) flux for the 3D hyperbolic diffusion equations.
+"""
 @inline function flux_godunov(u_ll, u_rr, orientation::Integer,
                               equations::HyperbolicDiffusionEquations3D)
     # Obtain left and right fluxes
diff --git a/src/equations/inviscid_burgers_1d.jl b/src/equations/inviscid_burgers_1d.jl
index afc6f9999c..2f9d3d1f6f 100644
--- a/src/equations/inviscid_burgers_1d.jl
+++ b/src/equations/inviscid_burgers_1d.jl
@@ -100,7 +100,14 @@ end
     return (abs(u[1]),)
 end
 
-# (Symmetric) Entropy Conserving flux
+@doc raw"""
+    flux_ec(u_ll, u_rr, orientation, equations::InviscidBurgersEquation1D)
+
+Entropy-conserving, symmetric flux for the inviscid Burgers' equation.
+```math
+F(u_L, u_R) = \frac{u_L^2 + u_L u_R + u_R^2}{6}
+```
+"""
 function flux_ec(u_ll, u_rr, orientation, equation::InviscidBurgersEquation1D)
     u_L = u_ll[1]
     u_R = u_rr[1]
@@ -108,8 +115,13 @@ function flux_ec(u_ll, u_rr, orientation, equation::InviscidBurgersEquation1D)
     return SVector((u_L^2 + u_L * u_R + u_R^2) / 6)
 end
 
-# See https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf ,
-# section 4.1.5 and especially equation (4.16).
+"""
+    flux_godunov(u_ll, u_rr, orientation, equations::InviscidBurgersEquation1D)
+
+Godunov (upwind) numerical flux for the inviscid Burgers' equation.
+See https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf ,
+section 4.1.5 and especially equation (4.16).
+"""
 function flux_godunov(u_ll, u_rr, orientation, equation::InviscidBurgersEquation1D)
     u_L = u_ll[1]
     u_R = u_rr[1]
diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl
index 920f862f43..b2323789dd 100644
--- a/src/equations/lattice_boltzmann_2d.jl
+++ b/src/equations/lattice_boltzmann_2d.jl
@@ -253,6 +253,12 @@ end
 #   λ_max =
 # end
 
+"""
+    flux_godunov(u_ll, u_rr, orientation, 
+                 equations::LatticeBoltzmannEquations2D)
+
+Godunov (upwind) flux for the 2D Lattice-Boltzmann equations.
+"""
 @inline function flux_godunov(u_ll, u_rr, orientation::Integer,
                               equations::LatticeBoltzmannEquations2D)
     if orientation == 1
diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl
index 1451eb7cbc..2211ac5bb7 100644
--- a/src/equations/lattice_boltzmann_3d.jl
+++ b/src/equations/lattice_boltzmann_3d.jl
@@ -243,6 +243,12 @@ end
 #   λ_max =
 # end
 
+"""
+    flux_godunov(u_ll, u_rr, orientation, 
+                 equations::LatticeBoltzmannEquations3D)
+
+Godunov (upwind) flux for the 3D Lattice-Boltzmann equations.
+"""
 @inline function flux_godunov(u_ll, u_rr, orientation::Integer,
                               equations::LatticeBoltzmannEquations3D)
     if orientation == 1 # x-direction
diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl
index 743d2df870..b1db57b0fd 100644
--- a/src/equations/linear_scalar_advection_1d.jl
+++ b/src/equations/linear_scalar_advection_1d.jl
@@ -141,8 +141,13 @@ end
     λ_max = abs(equation.advection_velocity[orientation])
 end
 
-# Essentially first order upwind, see e.g.
-# https://math.stackexchange.com/a/4355076/805029
+"""
+    flux_godunov(u_ll, u_rr, orientation, 
+                 equations::LinearScalarAdvectionEquation1D)
+
+Godunov (upwind) flux for the 1D linear scalar advection equation.
+Essentially first order upwind, see e.g. https://math.stackexchange.com/a/4355076/805029 .
+"""
 function flux_godunov(u_ll, u_rr, orientation::Int,
                       equation::LinearScalarAdvectionEquation1D)
     u_L = u_ll[1]
diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl
index 5e4f8463f5..3425870643 100644
--- a/src/equations/linear_scalar_advection_2d.jl
+++ b/src/equations/linear_scalar_advection_2d.jl
@@ -239,8 +239,13 @@ end
     return abs(a)
 end
 
-# Essentially first order upwind, see e.g.
-# https://math.stackexchange.com/a/4355076/805029
+"""
+    flux_godunov(u_ll, u_rr, orientation_or_normal_direction, 
+                 equations::LinearScalarAdvectionEquation2D)
+
+Godunov (upwind) flux for the 2D linear scalar advection equation.
+Essentially first order upwind, see e.g. https://math.stackexchange.com/a/4355076/805029 .
+"""
 function flux_godunov(u_ll, u_rr, orientation::Integer,
                       equation::LinearScalarAdvectionEquation2D)
     u_L = u_ll[1]
@@ -254,8 +259,6 @@ function flux_godunov(u_ll, u_rr, orientation::Integer,
     end
 end
 
-# Essentially first order upwind, see e.g.
-# https://math.stackexchange.com/a/4355076/805029
 function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
                       equation::LinearScalarAdvectionEquation2D)
     u_L = u_ll[1]
diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl
index 088f934cc3..689989d89c 100644
--- a/src/equations/linear_scalar_advection_3d.jl
+++ b/src/equations/linear_scalar_advection_3d.jl
@@ -158,8 +158,13 @@ end
     return abs(a)
 end
 
-# Essentially first order upwind, see e.g.
-# https://math.stackexchange.com/a/4355076/805029
+"""
+    flux_godunov(u_ll, u_rr, orientation_or_normal_direction, 
+                 equations::LinearScalarAdvectionEquation3D)
+
+Godunov (upwind) flux for the 3D linear scalar advection equation.
+Essentially first order upwind, see e.g. https://math.stackexchange.com/a/4355076/805029 .
+"""
 function flux_godunov(u_ll, u_rr, orientation::Integer,
                       equation::LinearScalarAdvectionEquation3D)
     u_L = u_ll[1]
@@ -173,8 +178,6 @@ function flux_godunov(u_ll, u_rr, orientation::Integer,
     end
 end
 
-# Essentially first order upwind, see e.g.
-# https://math.stackexchange.com/a/4355076/805029
 function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
                       equation::LinearScalarAdvectionEquation3D)
     u_L = u_ll[1]