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 OrdinaryDiffEqLowStorageRK
using 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-2
equations_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, n_cells_max = 30_000) # set maximum capacity of tree data structure
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: ………………………………… ViscousFormulationBassiRebay1 │
│ 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-6
sol = 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: ………………………………… ViscousFormulationBassiRebay1 │
│ 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 │
│ maximum #cells: ……………………………………… 30000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ 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.18,…iter=default_dt_factor_limiter) │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: ……………………………………………………… 1 │
│ threading backend: ……………………………… polyester │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
───────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 303ms / 98.1% 130KiB / 21.7%
Section ncalls time %tot avg alloc %tot avg
───────────────────────────────────────────────────────────────────────────────────────
parabolic rhs! 1.20k 222ms 74.6% 185μs 18.9KiB 66.9% 16.1B
calculate gradient 1.20k 101ms 34.1% 84.4μs 7.34KiB 26.0% 6.27B
volume integral 1.20k 64.1ms 21.6% 53.4μs 0.00B 0.0% 0.00B
surface integral 1.20k 10.3ms 3.5% 8.62μs 0.00B 0.0% 0.00B
Jacobian 1.20k 8.85ms 3.0% 7.38μs 0.00B 0.0% 0.00B
prolong2interfaces 1.20k 7.33ms 2.5% 6.11μs 0.00B 0.0% 0.00B
interface flux 1.20k 6.20ms 2.1% 5.17μs 0.00B 0.0% 0.00B
reset gradients 1.20k 1.42ms 0.5% 1.19μs 0.00B 0.0% 0.00B
~calculate gradient~ 1.20k 1.20ms 0.4% 1.00μs 7.34KiB 26.0% 6.27B
prolong2boundaries 1.20k 1.01ms 0.3% 840ns 0.00B 0.0% 0.00B
boundary flux 1.20k 743μs 0.3% 619ns 0.00B 0.0% 0.00B
prolong2mortars 1.20k 68.5μs 0.0% 57.1ns 0.00B 0.0% 0.00B
mortar flux 1.20k 52.5μs 0.0% 43.8ns 0.00B 0.0% 0.00B
volume integral 1.20k 52.3ms 17.6% 43.6μs 0.00B 0.0% 0.00B
calculate viscous fluxes 1.20k 38.3ms 12.9% 31.9μs 0.00B 0.0% 0.00B
prolong2interfaces 1.20k 7.33ms 2.5% 6.11μs 0.00B 0.0% 0.00B
transform variables 1.20k 6.60ms 2.2% 5.50μs 0.00B 0.0% 0.00B
interface flux 1.20k 6.24ms 2.1% 5.20μs 0.00B 0.0% 0.00B
surface integral 1.20k 4.78ms 1.6% 3.98μs 0.00B 0.0% 0.00B
~parabolic rhs!~ 1.20k 1.65ms 0.6% 1.37μs 11.5KiB 40.9% 9.84B
boundary flux 1.20k 978μs 0.3% 815ns 0.00B 0.0% 0.00B
prolong2boundaries 1.20k 772μs 0.3% 643ns 0.00B 0.0% 0.00B
reset ∂u/∂t 1.20k 684μs 0.2% 570ns 0.00B 0.0% 0.00B
Jacobian 1.20k 668μs 0.2% 557ns 0.00B 0.0% 0.00B
prolong2mortars 1.20k 62.8μs 0.0% 52.3ns 0.00B 0.0% 0.00B
mortar flux 1.20k 56.9μs 0.0% 47.4ns 0.00B 0.0% 0.00B
source terms parabolic 1.20k 35.9μs 0.0% 29.9ns 0.00B 0.0% 0.00B
rhs! 1.20k 75.4ms 25.4% 62.8μs 9.33KiB 33.1% 7.96B
volume integral 1.20k 48.2ms 16.2% 40.1μs 0.00B 0.0% 0.00B
interface flux 1.20k 13.7ms 4.6% 11.5μs 0.00B 0.0% 0.00B
surface integral 1.20k 4.81ms 1.6% 4.01μs 0.00B 0.0% 0.00B
prolong2interfaces 1.20k 3.86ms 1.3% 3.21μs 0.00B 0.0% 0.00B
boundary flux 1.20k 1.59ms 0.5% 1.33μs 0.00B 0.0% 0.00B
~rhs!~ 1.20k 1.28ms 0.4% 1.07μs 9.33KiB 33.1% 7.96B
Jacobian 1.20k 680μs 0.2% 567ns 0.00B 0.0% 0.00B
reset ∂u/∂t 1.20k 616μs 0.2% 513ns 0.00B 0.0% 0.00B
prolong2boundaries 1.20k 525μs 0.2% 437ns 0.00B 0.0% 0.00B
prolong2mortars 1.20k 59.0μs 0.0% 49.2ns 0.00B 0.0% 0.00B
mortar flux 1.20k 50.7μs 0.0% 42.3ns 0.00B 0.0% 0.00B
source terms 1.20k 35.9μs 0.0% 29.9ns 0.00B 0.0% 0.00B
───────────────────────────────────────────────────────────────────────────────────────We can now visualize the solution, which develops a boundary layer at the outflow boundaries.
using Plots
plot(sol)
Package versions
These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["Trixi", "OrdinaryDiffEqLowStorageRK", "Plots"],
mode = PKGMODE_MANIFEST)Julia Version 1.10.10
Commit 95f30e51f41 (2025-06-27 09:51 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 v1.10.0
[91a5bcdd] Plots v1.41.4
[a7f1ee26] Trixi v0.13.23 `~/work/Trixi.jl/Trixi.jl`This page was generated using Literate.jl.