trixi-framework · ranocha · Mar 26, 2026 · Mar 17, 2026 · Mar 17, 2026 · Mar 19, 2026
diff --git a/NEWS.md b/NEWS.md
@@ -8,7 +8,11 @@ for human readability.
 ## Changes in the v0.15 lifecycle
 
 #### Added
-
+- The word "viscous" is now used only where it refers specifically to fluid viscosity.
+The word "parabolic" is used in more general contexts.
+For example, `ViscousFormulationLocalDG()` is now `ParabolicFormulationLocalDG()`.
+In particular, viscosity is no longer used as a proxy for any parabolic/diffusive process such as heat conduction.
+For consistency, `cfl_advective` and `cfl_diffusive` have also been renamed `cfl_hyperbolic` and `cfl_parabolic` ([#2868]).
 - It is now possible to use `ViscousFormulationLocalDG()` as the `solver_parabolic` for non-conforming `P4estMesh`es.
 This is useful for (locally) diffusion-dominated problems.
 This enables in particular adaptive mesh refinement for that solver-mesh combination ([#2712]).

diff --git a/docs/literate/src/files/adding_new_parabolic_terms.jl b/docs/literate/src/files/adding_new_parabolic_terms.jl
@@ -35,14 +35,14 @@ function varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiff
     return varnames(variable_mapping, equations_parabolic.equations_hyperbolic)
 end
 
-# Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system
+# Next, we define the parabolic flux function. We assume that the mixed hyperbolic-parabolic system
 # is of the form
 # ```math
 # \partial_t u(t,x) + \partial_x (f_1(u) - g_1(u, \nabla u))
 #                   + \partial_y (f_2(u) - g_2(u, \nabla u)) = 0
 # ```
 # where ``f_1(u)``, ``f_2(u)`` are the hyperbolic fluxes and ``g_1(u, \nabla u)``, ``g_2(u, \nabla u)`` denote
-# the viscous fluxes. For anisotropic diffusion, the viscous fluxes are the first and second components
+# the parabolic fluxes. For anisotropic diffusion, the parabolic fluxes are the first and second components
 # of the matrix-vector product involving `diffusivity` and the gradient vector.
 #
 # Here, we specialize the flux to our new parabolic equation type `ConstantAnisotropicDiffusion2D`.
@@ -66,12 +66,12 @@ end
 # \begin{aligned}
 # \bm{q} &= \nabla u \\
 # \bm{\sigma} &= \begin{pmatrix} g_1(u, \bm{q}) \\ g_2(u, \bm{q}) \end{pmatrix} \\
-# \text{viscous contribution } &= \nabla \cdot \bm{\sigma}
+# \text{parabolic contribution} &= \nabla \cdot \bm{\sigma}
 # \end{aligned}
 # ```
 #
 # Boundary data must be specified for all spatial derivatives, e.g., for both the gradient
-# equation ``\bm{q} = \nabla u`` and the divergence of the viscous flux
+# equation ``\bm{q} = \nabla u`` and the divergence of the parabolic flux
 # ``\nabla \cdot \bm{\sigma}``. We account for this by introducing internal `Gradient`
 # and `Divergence` types which are used to dispatch on each type of boundary condition.
 #
@@ -98,7 +98,7 @@ end
     return boundary_condition.boundary_value
 end
 
-# While the gradient acts on the solution `u`, the divergence acts on the viscous flux ``\bm{\sigma}``.
+# While the gradient acts on the solution `u`, the divergence acts on the parabolic flux ``\bm{\sigma}``.
 # Thus, we have to supply boundary data for the `Divergence` operator that corresponds to ``\bm{\sigma}``.
 # However, we've already imposed boundary data on `u` for a Dirichlet boundary condition, and imposing
 # boundary data for ``\bm{\sigma}`` might overconstrain our problem.
@@ -119,7 +119,7 @@ end
 # ### A note on the choice of gradient variables
 #
 # It is often simpler to transform the solution variables (and solution gradients) to another set of
-# variables prior to computing the viscous fluxes (see [`CompressibleNavierStokesDiffusion2D`](@ref)
+# variables prior to computing the parabolic fluxes (see [`CompressibleNavierStokesDiffusion2D`](@ref)
 # for an example of this). If this is done, then the boundary condition for the `Gradient` operator
 # should be modified accordingly as well.
 #
@@ -182,7 +182,7 @@ plot(sol)
 # To be able to do so, we need to define [`max_diffusivity`](@ref) and 
 # [`have_constant_diffusivity`](@ref) for the new parabolic terms.
 # In Trixi.jl, currently only the standard Laplace Diffusion and Compressible Navier-Stokes-Fourier 
-# viscous terms are implemented.
+# parabolic terms are implemented.
 # Since these equations have **isotropic** diffusivity, i.e., direction-independent coefficients,
 # [`max_diffusivity`](@ref) is expected to return a scalar value.
 #
@@ -204,11 +204,11 @@ end
     return lambda_max()
 end
 
-# We supply now the advective(hyperbolic) and diffusive(parabolic) CFL numbers
-cfl_advective = 2.0 # Not restrictive for this example
-cfl_diffusive = 0.21 # Restricts the timestep
-stepsize_callback = StepsizeCallback(cfl = cfl_advective,
-                                     cfl_diffusive = cfl_diffusive)
+# We now supply the hyperbolic and parabolic CFL numbers
+cfl_hyperbolic = 2.0 # Not restrictive for this example
+cfl_parabolic = 0.21 # Restricts the timestep
+stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic,
+                                     cfl_parabolic = cfl_parabolic)
 
 # Add the stepsize callback to the existing callbacks
 callbacks = CallbackSet(SummaryCallback(), stepsize_callback);

diff --git a/docs/literate/src/files/parabolic_source_terms.jl b/docs/literate/src/files/parabolic_source_terms.jl
@@ -62,7 +62,7 @@ end
 # to OrdinaryDiffEq.jl. 
 # 
 # Note that for this problem, since viscosity `nu` is relatively large, we utilize 
-# `ViscousFormulationLocalDG` instead of the default `ViscousFormulationBassiRebay1` 
+# `ParabolicFormulationLocalDG` instead of the default `ParabolicFormulationBassiRebay1` 
 # parabolic solver, since the Bassi-Rebay 1 formulation is not accurate when the 
 # diffusivity is large relative to the mesh size. 
 
@@ -76,7 +76,7 @@ boundary_conditions_parabolic = boundary_condition_periodic
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationLocalDG(),
+                                             solver_parabolic = ParabolicFormulationLocalDG(),
                                              source_terms = source_terms,
                                              source_terms_parabolic = source_terms_parabolic,
                                              boundary_conditions = (boundary_conditions,
@@ -89,10 +89,10 @@ ode = semidiscretize(semi, tspan)
 # stable time-step is $O(h^2)$ due to the dominant parabolic term. We enforce this more stringent
 # parabolic CFL condition using a diffusion-aware `StepsizeCallback`. 
 
-cfl_advective = 0.5
-cfl_diffusive = 0.05
-stepsize_callback = StepsizeCallback(cfl = cfl_advective,
-                                     cfl_diffusive = cfl_diffusive)
+cfl_hyperbolic = 0.5
+cfl_parabolic = 0.05
+stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic,
+                                     cfl_parabolic = cfl_parabolic)
 callbacks = CallbackSet(SummaryCallback(), stepsize_callback)
 sol = solve(ode, RDPK3SpFSAL35(); adaptive = false, dt = stepsize_callback(ode),
             ode_default_options()..., callback = callbacks)
diff --git a/examples/dgmulti_1d/elixir_advection_diffusion_gradient_source_terms.jl b/examples/dgmulti_1d/elixir_advection_diffusion_gradient_source_terms.jl
@@ -58,10 +58,10 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype
 
 alive_callback = AliveCallback(analysis_interval = 100)
 
-cfl_advective = 0.5   # Not restrictive for this example
-cfl_diffusive = 0.025 # Restricts the timestep
-stepsize_callback = StepsizeCallback(cfl = cfl_advective,
-                                     cfl_diffusive = cfl_diffusive)
+cfl_hyperbolic = 0.5   # Not restrictive for this example
+cfl_parabolic = 0.025 # Restricts the timestep
+stepsize_callback = StepsizeCallback(cfl = cfl_hyperbolic,
+                                     cfl_parabolic = cfl_parabolic)
 
 callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
                         stepsize_callback)

diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/examples/dgmulti_2d/elixir_advection_diffusion.jl
@@ -34,7 +34,7 @@ boundary_conditions = (; left = boundary_condition_left,
                        top = boundary_condition_do_nothing,
                        right = boundary_condition_do_nothing)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; left = boundary_condition_left,
                                  bottom = boundary_condition_zero,
                                  top = boundary_condition_zero,

diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_nonperiodic.jl
@@ -53,7 +53,7 @@ boundary_conditions = (; left = boundary_condition,
                        bottom = boundary_condition,
                        right = boundary_condition_do_nothing)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; entire_boundary = boundary_condition)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),

diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl
@@ -198,7 +198,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to
 # define inviscid boundary conditions
 boundary_conditions = (; top_bottom = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),

diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl
@@ -206,7 +206,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to
 # define inviscid boundary conditions
 boundary_conditions = (; top_bottom = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),

diff --git a/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl b/examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl
@@ -51,7 +51,7 @@ boundary_condition_cavity = BoundaryConditionNavierStokesWall(velocity_bc_cavity
 boundary_conditions = (; top = boundary_condition_slip_wall,
                        rest_of_boundary = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; top = boundary_condition_lid,
                                  rest_of_boundary = boundary_condition_cavity)
 

diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl
@@ -241,7 +241,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to
 # define inviscid boundary conditions
 boundary_conditions = (; top_bottom = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),

diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl
@@ -249,7 +249,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to
 # define inviscid boundary conditions
 boundary_conditions = (; top_bottom = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; top_bottom = boundary_condition_top_bottom)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),

diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_nonperiodic_amr.jl
@@ -50,7 +50,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition)
 semi = SemidiscretizationHyperbolicParabolic(mesh,
                                              (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions_parabolic))
 

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
@@ -199,7 +199,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to
 boundary_conditions = (; y_neg = boundary_condition_slip_wall,
                        y_pos = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; y_neg = boundary_condition_top_bottom,
                                  y_pos = boundary_condition_top_bottom)
 

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
@@ -203,7 +203,7 @@ boundary_conditions = (; x_neg = boundary_condition_left_right,
                        y_neg = boundary_condition_slip_wall,
                        y_pos = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; x_neg = boundary_condition_left_right,
                                  x_pos = boundary_condition_left_right,
                                  y_neg = boundary_condition_top_bottom,

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_freestream_ldg.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_freestream_ldg.jl
@@ -19,7 +19,7 @@ initial_condition = initial_condition_const
 polydeg = 3
 solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs,
                volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
-solver_parabolic = ViscousFormulationLocalDG()
+solver_parabolic = ParabolicFormulationLocalDG()
 
 mu() = 0.5
 prandtl_number() = 0.72

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_shearlayer_nonconforming.jl
@@ -68,7 +68,7 @@ Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL)
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_condition_periodic,
                                                                     boundary_condition_periodic))
 

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_vortex_street.jl
@@ -119,8 +119,8 @@ function Trixi.get_node_variable(::Val{:vorticity}, u, mesh, equations, dg, cach
                             n_nodes, n_nodes, # equivalent: `ntuple(_ -> n_nodes, ndims(mesh))...,`
                             n_elements)
 
-    @unpack viscous_container = cache_parabolic
-    @unpack gradients = viscous_container
+    @unpack parabolic_container = cache_parabolic
+    @unpack gradients = parabolic_container
     gradients_x, gradients_y = gradients
 
     # We can accelerate the computation by thread-parallelizing the loop over elements

diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_amr_curved.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_amr_curved.jl
@@ -64,7 +64,7 @@ mesh = P4estMesh{3}(mesh_file, polydeg = 2,
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions))
 

diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonconforming.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonconforming.jl
@@ -70,7 +70,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition)
 semi = SemidiscretizationHyperbolicParabolic(mesh,
                                              (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions_parabolic))
 

diff --git a/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl b/examples/p4est_3d_dgsem/elixir_advection_diffusion_nonperiodic.jl
@@ -58,7 +58,7 @@ boundary_conditions = BoundaryConditionDirichlet(initial_condition)
 semi = SemidiscretizationHyperbolicParabolic(mesh,
                                              (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions))
 
@@ -76,7 +76,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
 stepsize_callback = StepsizeCallback(cfl = 1.6,
-                                     cfl_diffusive = 0.25)
+                                     cfl_parabolic = 0.25)
 
 callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
                         stepsize_callback)

diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl
@@ -67,7 +67,7 @@ mesh = P4estMesh(trees_per_dimension, polydeg = 3,
 
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_condition_periodic,
                                                                     boundary_condition_periodic))
 

diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl
@@ -243,7 +243,7 @@ boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_to
 boundary_conditions = (; y_neg = boundary_condition_slip_wall,
                        y_pos = boundary_condition_slip_wall)
 
-# define viscous boundary conditions
+# define parabolic boundary conditions
 boundary_conditions_parabolic = (; y_neg = boundary_condition_top_bottom,
                                  y_pos = boundary_condition_top_bottom)
 

diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_freestream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_freestream_boundaries.jl
@@ -19,7 +19,7 @@ initial_condition = initial_condition_const
 
 polydeg = 3
 solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
-solver_parabolic = ViscousFormulationBassiRebay1()
+solver_parabolic = ParabolicFormulationBassiRebay1()
 
 mu() = 0.5
 prandtl_number() = 0.72

diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_cfl.jl
@@ -62,10 +62,10 @@ analysis_callback = AnalysisCallback(semi, interval = 100)
 alive_callback = AliveCallback(analysis_interval = 100)
 
 # Stepsize callback which selects the timestep according to the most restrictive CFL condition.
-# For coarser grids, linear stability is governed by the advective CFL condition,
+# For coarser grids, linear stability is governed by the hyperbolic CFL condition,
 # while for high refinements the flow becomes diffusion-dominated.
 stepsize_callback = StepsizeCallback(cfl = 1.6,
-                                     cfl_diffusive = 0.3)
+                                     cfl_parabolic = 0.3)
 
 callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
                         stepsize_callback)

diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_dirichlet_amr.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_dirichlet_amr.jl
@@ -42,7 +42,7 @@ boundary_conditions_parabolic = BoundaryConditionDirichlet(initial_condition)
 semi = SemidiscretizationHyperbolicParabolic(mesh,
                                              (equations, equations_parabolic),
                                              initial_condition, solver;
-                                             solver_parabolic = ViscousFormulationBassiRebay1(),
+                                             solver_parabolic = ParabolicFormulationBassiRebay1(),
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions_parabolic))