13: Parabolic terms
Experimental support for parabolic diffusion terms is available in Trixi.jl. This demo illustrates parabolic terms for the advection-diffusion equation.
using OrdinaryDiffEqLowStorageRKusing TrixiSplitting a system into hyperbolic and parabolic parts
For a mixed hyperbolic-parabolic system, we represent the hyperbolic and parabolic parts of the system separately. We first define the hyperbolic (advection) part of the advection-diffusion equation.
advection_velocity = (1.5, 1.0)equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity);Next, we define the parabolic diffusion term. The constructor requires knowledge of equations_hyperbolic to be passed in because the LaplaceDiffusion2D applies diffusion to every variable of the hyperbolic system.
diffusivity = 5.0e-2equations_parabolic = LaplaceDiffusion2D(diffusivity, equations_hyperbolic);Boundary conditions
As with the equations, we define boundary conditions separately for the hyperbolic and parabolic part of the system. For this example, we impose inflow BCs for the hyperbolic system (no condition is imposed on the outflow), and we impose Dirichlet boundary conditions for the parabolic equations. Both BoundaryConditionDirichlet and BoundaryConditionNeumann are defined for LaplaceDiffusion2D.
The hyperbolic and parabolic boundary conditions are assumed to be consistent with each other.
boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])), y_neg = boundary_condition_zero_dirichlet, y_pos = boundary_condition_do_nothing, x_pos = boundary_condition_do_nothing)boundary_conditions_parabolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])), y_neg = boundary_condition_zero_dirichlet, y_pos = boundary_condition_zero_dirichlet, x_pos = boundary_condition_zero_dirichlet);Defining the solver and mesh
The process of creating the DG solver and mesh is the same as for a purely hyperbolic system of equations.
solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 4, periodicity = false)initial_condition = (x, t, equations) -> SVector(0.0);Semidiscretizing and solving
To semidiscretize a hyperbolic-parabolic system, we create a SemidiscretizationHyperbolicParabolic. This differs from a SemidiscretizationHyperbolic in that we pass in a Tuple containing both the hyperbolic and parabolic equation, as well as a Tuple containing the hyperbolic and parabolic boundary conditions.
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations_hyperbolic, equations_parabolic), initial_condition, solver; boundary_conditions = (boundary_conditions_hyperbolic, boundary_conditions_parabolic))┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolicParabolic │
│ ═════════════════════════════════════ │
│ #spatial dimensions: ……………………………………… 2 │
│ mesh: ……………………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2, Float64}} with length 341 │
│ hyperbolic equations: …………………………………… LinearScalarAdvectionEquation2D │
│ parabolic equations: ……………………………………… LaplaceDiffusion2D │
│ initial condition: …………………………………………… #7 │
│ source terms: ………………………………………………………… nothing │
│ source terms parabolic: ……………………………… nothing │
│ solver: ………………………………………………………………………… DG │
│ parabolic solver: ……………………………………………… ParabolicFormulationBassiRebay1 │
│ total #DOFs per field: ………………………………… 4096 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘The rest of the code is identical to the hyperbolic case. We create a system of ODEs through semidiscretize, defining callbacks, and then passing the system to OrdinaryDiffEq.jl.
tspan = (0.0, 1.5)ode = semidiscretize(semi, tspan)callbacks = CallbackSet(SummaryCallback())time_int_tol = 1.0e-6sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, ode_default_options()..., callback = callbacks);
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolicParabolic │
│ ═════════════════════════════════════ │
│ #spatial dimensions: ……………………………………… 2 │
│ mesh: ……………………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2, Float64}} with length 341 │
│ hyperbolic equations: …………………………………… LinearScalarAdvectionEquation2D │
│ parabolic equations: ……………………………………… LaplaceDiffusion2D │
│ initial condition: …………………………………………… #7 │
│ source terms: ………………………………………………………… nothing │
│ source terms parabolic: ……………………………… nothing │
│ solver: ………………………………………………………………………… DG │
│ parabolic solver: ……………………………………………… ParabolicFormulationBassiRebay1 │
│ total #DOFs per field: ………………………………… 4096 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2, Float64}} │
│ ═════════════════════════════════════════ │
│ center: ………………………………………………………………………… [0.0, 0.0] │
│ length: ………………………………………………………………………… 2.0 │
│ periodicity: …………………………………………………………… (false, false) │
│ current #cells: …………………………………………………… 341 │
│ #leaf-cells: …………………………………………………………… 256 │
│ current capacity: ……………………………………………… 341 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation2D │
│ ═══════════════════════════════ │
│ #variables: ……………………………………………………………… 1 │
│ │ variable 1: ………………………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: …………………………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: ………………………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ……………………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: …………………………………………………… FluxLaxFriedrichs(max_abs_speed) │
│ volume integral: ………………………………………………… VolumeIntegralWeakForm │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ……………………………………………………………… 0.0 │
│ Final time: ……………………………………………………………… 1.5 │
│ time integrator: ………………………………………………… RDPK3SpFSAL49 │
│ adaptive: …………………………………………………………………… true │
│ abstol: ………………………………………………………………………… 1.0e-6 │
│ reltol: ………………………………………………………………………… 1.0e-6 │
│ controller: ……………………………………………………………… PIDController(beta=(0.38, -0.…er=default_dt_factor_limiter) │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: …………………………………………………………………… 1 │
│ threading backend: …………………………………………… polyester │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
Trixi.jl
───────────────────────────────────────────────────────────────────────────────────────────
Time Allocations
────────────────────── ────────────────────────
Tot / % measured: 101ms / 94.4% 115KiB / 9.8%
───────────────────────────────────── ────────────────────── ────────────────────────
Section ncalls time %tot avg alloc %tot avg
───────────────────────────────────────────────────────────────────────────────────────────
parabolic rhs! 1.20k 66.3ms 69.7% 55.3μs 7.19KiB 64.0% 6.13B
├─ calculate gradient 1.20k 31.9ms 33.6% 26.6μs 2.62KiB 23.4% 2.24B
│ ├─ volume integral 1.20k 20.1ms 21.1% 16.7μs ∅ ∅ ∅
│ ├─ surface integral 1.20k 2.67ms 2.8% 2.23μs ∅ ∅ ∅
│ ├─ interface flux 1.20k 2.23ms 2.3% 1.86μs ∅ ∅ ∅
│ ├─ prolong2interfaces 1.20k 2.03ms 2.1% 1.69μs ∅ ∅ ∅
│ ├─ reset gradients 1.20k 1.91ms 2.0% 1.59μs ∅ ∅ ∅
│ ├─ ~calculate gradient~ 1.20k 1.26ms 1.3% 1.05μs 2.62KiB 23.4% 2.24B
│ ├─ Jacobian 1.20k 999μs 1.1% 833ns ∅ ∅ ∅
│ ├─ prolong2boundaries 1.20k 371μs 0.4% 309ns ∅ ∅ ∅
│ ├─ boundary flux 1.20k 213μs 0.2% 178ns ∅ ∅ ∅
│ ├─ prolong2mortars 1.20k 88.7μs 0.1% 73.9ns ∅ ∅ ∅
│ └─ mortar flux 1.20k 76.3μs 0.1% 63.5ns ∅ ∅ ∅
├─ volume integral 1.20k 18.5ms 19.4% 15.4μs ∅ ∅ ∅
├─ calculate parabolic fluxes 1.20k 3.53ms 3.7% 2.94μs ∅ ∅ ∅
├─ surface integral 1.20k 2.44ms 2.6% 2.03μs ∅ ∅ ∅
├─ interface flux 1.20k 2.20ms 2.3% 1.83μs ∅ ∅ ∅
├─ prolong2interfaces 1.20k 2.05ms 2.2% 1.71μs ∅ ∅ ∅
├─ ~parabolic rhs!~ 1.20k 1.86ms 2.0% 1.55μs 4.56KiB 40.6% 3.89B
├─ transform variables 1.20k 1.64ms 1.7% 1.37μs ∅ ∅ ∅
├─ reset ∂u/∂t 1.20k 633μs 0.7% 527ns ∅ ∅ ∅
├─ boundary flux 1.20k 519μs 0.5% 433ns ∅ ∅ ∅
├─ Jacobian 1.20k 491μs 0.5% 409ns ∅ ∅ ∅
├─ prolong2boundaries 1.20k 371μs 0.4% 309ns ∅ ∅ ∅
├─ prolong2mortars 1.20k 75.4μs 0.1% 62.8ns ∅ ∅ ∅
├─ mortar flux 1.20k 64.9μs 0.1% 54.1ns ∅ ∅ ∅
└─ source terms parabolic 1.20k 34.8μs 0.0% 29.0ns ∅ ∅ ∅
rhs_hyperbolic! 1.20k 28.8ms 30.3% 24.0μs 4.05KiB 36.0% 3.45B
├─ volume integral 1.20k 17.6ms 18.5% 14.7μs ∅ ∅ ∅
├─ interface flux 1.20k 3.02ms 3.2% 2.52μs ∅ ∅ ∅
├─ surface integral 1.20k 2.46ms 2.6% 2.05μs ∅ ∅ ∅
├─ prolong2interfaces 1.20k 2.17ms 2.3% 1.81μs ∅ ∅ ∅
├─ ~rhs_hyperbolic!~ 1.20k 1.43ms 1.5% 1.19μs 4.05KiB 36.0% 3.45B
├─ boundary flux 1.20k 642μs 0.7% 535ns ∅ ∅ ∅
├─ reset ∂u/∂t 1.20k 561μs 0.6% 468ns ∅ ∅ ∅
├─ Jacobian 1.20k 490μs 0.5% 408ns ∅ ∅ ∅
├─ prolong2boundaries 1.20k 245μs 0.3% 204ns ∅ ∅ ∅
├─ prolong2mortars 1.20k 70.0μs 0.1% 58.4ns ∅ ∅ ∅
├─ mortar flux 1.20k 61.9μs 0.1% 51.6ns ∅ ∅ ∅
└─ source terms 1.20k 35.0μs 0.0% 29.2ns ∅ ∅ ∅
───────────────────────────────────────────────────────────────────────────────────────────We can now visualize the solution, which develops a boundary layer at the outflow boundaries.
using Plotsplot(sol)
Package versions
These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.status(["Trixi", "OrdinaryDiffEqLowStorageRK", "Plots"], mode = PKGMODE_MANIFEST)Julia Version 1.10.12
Commit d93beab124c (2026-08-15 10:29 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 7763 64-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/Trixi.jl/Trixi.jl/docs/Manifest.toml`
[b0944070] OrdinaryDiffEqLowStorageRK v3.4.0
[91a5bcdd] Plots v1.41.7
[a7f1ee26] Trixi v0.17.15-DEV `~/work/Trixi.jl/Trixi.jl`This page was generated using Literate.jl.