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 Trixi

Splitting 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                                                          │
│ 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                                                          │
│ 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_naive)                           │
│ 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                                                                │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
───────────────────────────────────────────────────────────────────────────────────────
              Trixi.jl                        Time                    Allocations
                                     ───────────────────────   ────────────────────────
          Tot / % measured:               295ms /  98.0%            129KiB /  21.3%

Section                      ncalls     time    %tot     avg     alloc    %tot      avg
───────────────────────────────────────────────────────────────────────────────────────
parabolic rhs!                1.20k    214ms   74.0%   178μs   18.1KiB   66.0%    15.5B
  calculate gradient          1.20k    100ms   34.4%  82.9μs   7.34KiB   26.7%    6.27B
    volume integral           1.20k   61.2ms   21.2%  51.0μs     0.00B    0.0%    0.00B
    surface integral          1.20k   10.6ms    3.7%  8.80μs     0.00B    0.0%    0.00B
    Jacobian                  1.20k   8.91ms    3.1%  7.42μs     0.00B    0.0%    0.00B
    prolong2interfaces        1.20k   7.34ms    2.5%  6.12μs     0.00B    0.0%    0.00B
    interface flux            1.20k   6.90ms    2.4%  5.75μs     0.00B    0.0%    0.00B
    reset gradients           1.20k   1.54ms    0.5%  1.28μs     0.00B    0.0%    0.00B
    ~calculate gradient~      1.20k   1.15ms    0.4%   957ns   7.34KiB   26.7%    6.27B
    prolong2boundaries        1.20k    997μs    0.3%   831ns     0.00B    0.0%    0.00B
    boundary flux             1.20k    765μs    0.3%   638ns     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   52.0μs    0.0%  43.4ns     0.00B    0.0%    0.00B
  volume integral             1.20k   48.3ms   16.7%  40.2μs     0.00B    0.0%    0.00B
  calculate viscous fluxes    1.20k   36.1ms   12.5%  30.1μs     0.00B    0.0%    0.00B
  prolong2interfaces          1.20k   7.62ms    2.6%  6.35μs     0.00B    0.0%    0.00B
  transform variables         1.20k   6.52ms    2.3%  5.43μs     0.00B    0.0%    0.00B
  interface flux              1.20k   6.30ms    2.2%  5.25μs     0.00B    0.0%    0.00B
  surface integral            1.20k   4.74ms    1.6%  3.95μs     0.00B    0.0%    0.00B
  ~parabolic rhs!~            1.20k   1.47ms    0.5%  1.23μs   10.8KiB   39.3%    9.21B
  boundary flux               1.20k    993μs    0.3%   828ns     0.00B    0.0%    0.00B
  prolong2boundaries          1.20k    753μs    0.3%   627ns     0.00B    0.0%    0.00B
  reset ∂u/∂t                 1.20k    695μs    0.2%   579ns     0.00B    0.0%    0.00B
  Jacobian                    1.20k    655μs    0.2%   546ns     0.00B    0.0%    0.00B
  prolong2mortars             1.20k   63.5μs    0.0%  52.9ns     0.00B    0.0%    0.00B
  mortar flux                 1.20k   46.2μs    0.0%  38.5ns     0.00B    0.0%    0.00B
rhs!                          1.20k   75.2ms   26.0%  62.7μs   9.33KiB   34.0%    7.96B
  volume integral             1.20k   48.1ms   16.6%  40.1μs     0.00B    0.0%    0.00B
  interface flux              1.20k   13.7ms    4.8%  11.4μs     0.00B    0.0%    0.00B
  surface integral            1.20k   4.75ms    1.6%  3.96μs     0.00B    0.0%    0.00B
  prolong2interfaces          1.20k   3.79ms    1.3%  3.15μs     0.00B    0.0%    0.00B
  boundary flux               1.20k   1.56ms    0.5%  1.30μs     0.00B    0.0%    0.00B
  ~rhs!~                      1.20k   1.31ms    0.5%  1.09μs   9.33KiB   34.0%    7.96B
  Jacobian                    1.20k    706μs    0.2%   589ns     0.00B    0.0%    0.00B
  reset ∂u/∂t                 1.20k    580μs    0.2%   484ns     0.00B    0.0%    0.00B
  prolong2boundaries          1.20k    538μs    0.2%   449ns     0.00B    0.0%    0.00B
  prolong2mortars             1.20k   58.8μs    0.0%  49.0ns     0.00B    0.0%    0.00B
  mortar flux                 1.20k   47.1μs    0.0%  39.2ns     0.00B    0.0%    0.00B
  source terms                1.20k   36.5μs    0.0%  30.4ns     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)
Example block output

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.8
Commit 4c16ff44be8 (2025-01-22 10:06 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.2.1
  [91a5bcdd] Plots v1.40.9
  [a7f1ee26] Trixi v0.11.1-DEV `~/work/Trixi.jl/Trixi.jl`

This page was generated using Literate.jl.