(Linear) Matrix operator for efficient IMEX

Project.toml

@@ -70,14 +70,14 @@ TrixiSparseConnectivityTracerExt = "SparseConnectivityTracer"
 
 [compat]
 Accessors = "0.1.36"
-Adapt = "4"
-CUDA = "5.8"
+Adapt = "4.1"
+CUDA = "5.8.2"
 CodeTracking = "1.0.5, 2, 3"
 ConstructionBase = "1.5"
 Convex = "0.16"
 DataStructures = "0.18.15, 0.19"
 DelimitedFiles = "1"
-DiffEqBase = "6.155.2"
+DiffEqBase = "6.174"
 DiffEqCallbacks = "2.35, 3, 4"
 Downloads = "1.6"
 ECOS = "1.1.2"
@@ -89,7 +89,7 @@ KernelAbstractions = "0.9.36"
 LinearAlgebra = "1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.171"
-MPI = "0.20.6"
+MPI = "0.20.22"
 Makie = "0.21, 0.22, 0.23, 0.24"
 MuladdMacro = "0.2.4"
 NLsolve = "4.5.1"
@@ -104,7 +104,7 @@ RecipesBase = "1.3.4"
 RecursiveArrayTools = "3.31.1"
 Reexport = "1.2"
 Requires = "1.3"
-SciMLBase = "2.67.0"
+SciMLBase = "2.92.0"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 SparseConnectivityTracer = "1.0.1"

examples/tree_2d_dgsem/elixir_advection_diffusion_imex_operator.jl

@@ -0,0 +1,74 @@
+using OrdinaryDiffEqSDIRK
+using SciMLBase, SciMLOperators, SparseArrays
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+advection_velocity = (1.5, 1.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+diffusivity() = 0.5f0
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                periodicity = true,
+                n_cells_max = 30_000)
+
+function initial_condition_diffusive_convergence_test(x, t,
+                                                      equation::LinearScalarAdvectionEquation2D)
+    # Store translated coordinate for easy use of exact solution
+    RealT = eltype(x)
+    x_trans = x - equation.advection_velocity * t
+
+    nu = diffusivity()
+    c = 1
+    A = 0.5f0
+    L = 2
+    f = 1.0f0 / L
+    omega = 2 * convert(RealT, pi) * f
+    scalar = c + A * sin(omega * sum(x_trans)) * exp(-2 * nu * omega^2 * t)
+    return SVector(scalar)
+end
+initial_condition = initial_condition_diffusive_convergence_test
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE setup
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan) # For accessing the initial condition u0 only
+
+D_map, _ = linear_structure_parabolic(semi)
+# Cannot directly construct `MatrixOperator` from `LinearMap`, need detour via sparse matrix
+D_op = MatrixOperator(sparse(D_map))
+
+split_func = SplitFunction(D_op, Trixi.rhs!)
+ode_operator = SplitODEProblem{true}(split_func, ode.u0, tspan, semi)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode_operator, KenCarp4();
+            adaptive = false, dt = 1e-2,
+            ode_default_options()..., callback = callbacks)

src/Trixi.jl

@@ -253,6 +253,7 @@ export entropy, energy_total, energy_kinetic, energy_internal,
        enstrophy, vorticity
 export lake_at_rest_error
 export ncomponents, eachcomponent
+export have_constant_speed
 
 export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
        P4estMeshView, P4estMeshCubedSphere, T8codeMesh
@@ -321,7 +322,7 @@ export ode_norm, ode_unstable_check
 
 export convergence_test,
        jacobian_fd, jacobian_ad_forward, jacobian_ad_forward_parabolic,
-       linear_structure
+       linear_structure, linear_structure_parabolic
 
 export DGMulti, DGMultiBasis, estimate_dt, DGMultiMesh, GaussSBP
 

src/equations/acoustic_perturbation_2d.jl

@@ -399,6 +399,24 @@ end
     return SVector(diss[1], diss[2], diss[3], z, z, z, z)
 end
 
+"""
+    have_constant_speed(::AcousticPerturbationEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref).
+
+The acoustic perturbation equations are in principle linear for **constant** mean flow fields.
+However, the mean flow variables are part of the solution vector in
+[`AcousticPerturbationEquations2D`](@ref) and only the **global** mean flow variables are constant,
+similar to the [`LinearizedEulerEquations2D`](@ref).
+
+Moreover, when coupling to the [`CompressibleEulerEquations2D`](@ref) equations via 
+[`SemidiscretizationEulerAcoustics`](@ref), the mean field variables are updated
+on the fly, see [`EulerAcousticsCouplingCallback`](@ref).
+
+# Returns
+- `False()`
+"""
 @inline have_constant_speed(::AcousticPerturbationEquations2D) = False()
 
 @inline function max_abs_speeds(u, equations::AcousticPerturbationEquations2D)

src/equations/compressible_navier_stokes.jl

@@ -110,13 +110,16 @@ dynamic_viscosity(u, mu::Real, equations) = mu
 dynamic_viscosity(u, mu::T, equations) where {T} = mu(u, equations)
 
 """
-    max_diffusivity(::AbstractCompressibleNavierStokesDiffusion)
+    have_constant_diffusivity(::AbstractCompressibleNavierStokesDiffusion)
 
 # Returns
 - `False()`
 
 Used in diffusive CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the
 diffusivity is not constant in space and that [`max_diffusivity`](@ref) needs to be computed
 at every node in every element.
+
+Also employed in [`linear_structure`](@ref) and [`linear_structure_parabolic`](@ref) to check
+if the diffusion term is linear in the variables/constant.
 """
 @inline have_constant_diffusivity(::AbstractCompressibleNavierStokesDiffusion) = False()

src/equations/equations.jl

@@ -362,6 +362,17 @@ See also the test section P4estMesh2D with combine_conservative_and_nonconservat
 """
 combine_conservative_and_nonconservative_fluxes(flux, ::AbstractEquations) = False()
 
+"""
+    have_constant_speed(::AbstractEquations)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref).
+
+This is the default fallback for nonlinear equations.
+
+# Returns
+- `False()`
+"""
 have_constant_speed(::AbstractEquations) = False()
 
 """

src/equations/equations_parabolic.jl

@@ -12,14 +12,17 @@ abstract type AbstractLaplaceDiffusion{NDIMS, NVARS} <:
               AbstractEquationsParabolic{NDIMS, NVARS, GradientVariablesConservative} end
 
 """
-    max_diffusivity(::AbstractLaplaceDiffusion)
+    have_constant_diffusivity(::AbstractLaplaceDiffusion)
 
 # Returns
 - `True()`
 
 Used in diffusive CFL condition computation (see [`StepsizeCallback`](@ref)) to indicate that the
 diffusivity is constant in space and that [`max_diffusivity`](@ref) needs **not** to be re-computed
 at every node in every element.
+
+Also employed in [`linear_structure`](@ref) and [`linear_structure_parabolic`](@ref) to check
+if the diffusion term is linear in the variables/constant.
 """
 @inline have_constant_diffusivity(::AbstractLaplaceDiffusion) = True()
 

src/equations/hyperbolic_diffusion_1d.jl

@@ -172,6 +172,15 @@ end
     return sqrt(equations.nu * equations.inv_Tr)
 end
 
+"""
+    have_constant_speed(::HyperbolicDiffusionEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::HyperbolicDiffusionEquations1D) = True()
 
 @inline function max_abs_speeds(eq::HyperbolicDiffusionEquations1D)

src/equations/hyperbolic_diffusion_2d.jl

@@ -224,6 +224,15 @@ end
     return SVector(f1, f2, f3)
 end
 
+"""
+    have_constant_speed(::HyperbolicDiffusionEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::HyperbolicDiffusionEquations2D) = True()
 
 @inline function max_abs_speeds(eq::HyperbolicDiffusionEquations2D)

src/equations/hyperbolic_diffusion_3d.jl

@@ -223,6 +223,15 @@ Godunov (upwind) flux for the 3D hyperbolic diffusion equations.
     return SVector(f1, f2, f3, f4)
 end
 
+"""
+    have_constant_speed(::HyperbolicDiffusionEquations3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::HyperbolicDiffusionEquations3D) = True()
 
 @inline function max_abs_speeds(eq::HyperbolicDiffusionEquations3D)

src/equations/lattice_boltzmann_2d.jl

@@ -371,6 +371,15 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
     return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
 end
 
+"""
+    have_constant_speed(::LatticeBoltzmannEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LatticeBoltzmannEquations2D) = True()
 
 @inline function max_abs_speeds(equations::LatticeBoltzmannEquations2D)

src/equations/lattice_boltzmann_3d.jl

@@ -384,6 +384,15 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
     return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
 end
 
+"""
+    have_constant_speed(::LatticeBoltzmannEquations3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LatticeBoltzmannEquations3D) = True()
 
 @inline function max_abs_speeds(equations::LatticeBoltzmannEquations3D)

src/equations/linear_elasticity_1d.jl

@@ -97,6 +97,15 @@ end
     return SVector(f1, f2)
 end
 
+"""
+    have_constant_speed(::LinearElasticityEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearElasticityEquations1D) = True()
 
 @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,

src/equations/linear_scalar_advection_1d.jl

@@ -170,6 +170,15 @@ function flux_engquist_osher(u_ll, u_rr, orientation::Int,
                     abs(equation.advection_velocity[orientation]) * (u_R - u_L)))
 end
 
+"""
+    have_constant_speed(::LinearScalarAdvectionEquation1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearScalarAdvectionEquation1D) = True()
 
 @inline function max_abs_speeds(equation::LinearScalarAdvectionEquation1D)

src/equations/linear_scalar_advection_2d.jl

@@ -269,6 +269,15 @@ function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
     end
 end
 
+"""
+    have_constant_speed(::LinearScalarAdvectionEquation2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearScalarAdvectionEquation2D) = True()
 
 @inline function max_abs_speeds(equation::LinearScalarAdvectionEquation2D)

src/equations/linear_scalar_advection_3d.jl

@@ -188,6 +188,15 @@ function flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
     end
 end
 
+"""
+    have_constant_speed(::LinearScalarAdvectionEquation3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearScalarAdvectionEquation3D) = True()
 
 @inline function max_abs_speeds(equation::LinearScalarAdvectionEquation3D)

src/equations/linearized_euler_1d.jl

@@ -108,6 +108,15 @@ end
     return SVector(f1, f2, f3)
 end
 
+"""
+    have_constant_speed(::LinearizedEulerEquations1D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearizedEulerEquations1D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations1D)

src/equations/linearized_euler_2d.jl

@@ -144,6 +144,15 @@ end
     return SVector(f1, f2, f3, f4)
 end
 
+"""
+    have_constant_speed(::LinearizedEulerEquations2D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearizedEulerEquations2D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations2D)

src/equations/linearized_euler_3d.jl

@@ -180,6 +180,15 @@ end
     return SVector(f1, f2, f3, f4, f5)
 end
 
+"""
+    have_constant_speed(::LinearizedEulerEquations3D)
+
+Indicates whether the characteristic speeds are constant, i.e., independent of the solution.
+Queried in the timestep computation [`StepsizeCallback`](@ref) and [`linear_structure`](@ref).
+
+# Returns
+- `True()`
+"""
 @inline have_constant_speed(::LinearizedEulerEquations3D) = True()
 
 @inline function max_abs_speeds(equations::LinearizedEulerEquations3D)