(Linear) Matrix operator for efficient IMEX

Project.toml

@@ -69,58 +69,58 @@ TrixiNLsolveExt = "NLsolve"
 TrixiSparseConnectivityTracerExt = "SparseConnectivityTracer"
 
 [compat]
-Accessors = "0.1.36"
-Adapt = "4"
-CUDA = "5.8"
-CodeTracking = "1.0.5, 2, 3"
-ConstructionBase = "1.5"
+Accessors = "0.1.42"
+Adapt = "4.4"
+CUDA = "5.9"
+CodeTracking = "2.0.2"
+ConstructionBase = "1.6"
 Convex = "0.16"
-DataStructures = "0.18.15, 0.19"
-DelimitedFiles = "1"
-DiffEqBase = "6.155.2"
-DiffEqCallbacks = "2.35, 3, 4"
+DataStructures = "0.19.3"
+DelimitedFiles = "1.9.1"
+DiffEqBase = "6.192"
+DiffEqCallbacks = "4.10.1"
 Downloads = "1.6"
 ECOS = "1.1.2"
-EllipsisNotation = "1.0"
-FillArrays = "1.9"
-ForwardDiff = "0.10.36, 1"
-HDF5 = "0.16.10, 0.17"
-KernelAbstractions = "0.9.36"
+EllipsisNotation = "1.8"
+FillArrays = "1.15"
+ForwardDiff = "1.3"
+HDF5 = "0.17.2"
+KernelAbstractions = "0.9.39"
 LinearAlgebra = "1"
-LinearMaps = "2.7, 3.0"
-LoopVectorization = "0.12.171"
-MPI = "0.20.6"
-Makie = "0.21, 0.22, 0.23, 0.24"
+LinearMaps = "3.11.4"
+LoopVectorization = "0.12.173"
+MPI = "0.20.23"
+Makie = "0.24.7"
 MuladdMacro = "0.2.4"
 NLsolve = "4.5.1"
-Octavian = "0.3.28"
-OffsetArrays = "1.13"
-P4est = "0.4.12"
-Polyester = "=0.7.16, 0.7.18"
-PrecompileTools = "1.2"
-Preferences = "1.4"
+Octavian = "0.3.29"
+OffsetArrays = "1.17"
+P4est = "0.4.13"
+Polyester = "0.7.18"
+PrecompileTools = "1.2.1"
+Preferences = "1.5"
 Printf = "1"
 RecipesBase = "1.3.4"
-RecursiveArrayTools = "3.31.1"
-Reexport = "1.2"
-Requires = "1.3"
-SciMLBase = "2.67.0"
+RecursiveArrayTools = "3.39"
+Reexport = "1.2.2"
+Requires = "1.3.1"
+SciMLBase = "2.127"
 SimpleUnPack = "1.1"
-SparseArrays = "1"
+SparseArrays = "1.10"
 SparseConnectivityTracer = "1.0.1"
-StableRNGs = "1.0.2"
-StartUpDG = "1.1.5"
-Static = "1.1.1"
-StaticArrayInterface = "1.5.1"
-StaticArrays = "1.9"
+StableRNGs = "1.0.4"
+StartUpDG = "1.3.3"
+Static = "1.3.1"
+StaticArrayInterface = "1.8"
+StaticArrays = "1.9.15"
 StrideArrays = "0.1.29"
-StructArrays = "0.6.20, 0.7"
-SummationByPartsOperators = "0.5.52"
+StructArrays = "0.7.2"
+SummationByPartsOperators = "0.5.86"
 T8code = "0.7.4"
-TimerOutputs = "0.5.25"
-Triangulate = "2.2, 3"
+TimerOutputs = "0.5.29"
+Triangulate = "3"
 TriplotBase = "0.1"
-TriplotRecipes = "0.1"
-TrixiBase = "0.1.6"
+TriplotRecipes = "0.1.2"
+TrixiBase = "0.1.8"
 UUIDs = "1.6"
 julia = "1.10"

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() = 5.0e-1
+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

@@ -320,7 +320,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/semidiscretization/semidiscretization.jl

@@ -245,7 +245,8 @@ and a vector `b`:
 ```math
 \\partial_t u(t) = A u(t) - b.
 ```
-Works only for linear equations, i.e., equations with `have_constant_speed(equations) == True()`.
+Works only for linear equations, i.e.,
+equations which `have_constant_speed(equations) == True()`.
 
 This has the benefit of greatly reduced memory consumption compared to constructing
 the full system matrix explicitly, as done for instance in
@@ -268,6 +269,7 @@ function linear_structure(semi::AbstractSemidiscretization;
     u_ode .= zero(eltype(u_ode))
     rhs!(du_ode, u_ode, semi, t0)
     b = -du_ode
+
     # Create a copy of `b` used internally to extract the linear part of `semi`.
     # This is necessary to get everything correct when the user updates the
     # returned vector `b`.

src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl

@@ -394,7 +394,8 @@ and a vector `b`:
 ```math
 \\partial_t u(t) = A u(t) - b.
 ```
-Works only for linear equations, i.e., equations with `have_constant_speed(equations) == True()`.
+Works only for linear equations, i.e.,
+equations which `have_constant_speed(equations) == True()`.
 
 This has the benefit of greatly reduced memory consumption compared to constructing
 the full system matrix explicitly, as done for instance in
@@ -457,6 +458,70 @@ function _jacobian_ad_forward(semi::SemidiscretizationHyperbolicParabolic, t0, u
     return J
 end
 
+"""
+    linear_structure_parabolic(semi::SemidiscretizationHyperbolicParabolic;
+                               t0 = zero(real(semi)))
+
+Wraps the **parabolic part** right-hand side operator of the hyperbolic-parabolic semidiscretization `semi`
+at time `t0` as an affine-linear operator given by a linear operator `A`
+and a vector `b`:
+```math
+\\partial_t u(t) = A u(t) - b.
+```
+Works only for linear parabolic equations, i.e.,
+equations which `have_constant_diffusivity(equations_parabolic) == True()`.
+
+This has the benefit of greatly reduced memory consumption compared to constructing
+the full system matrix explicitly, as done for instance in
+[`jacobian_ad_forward_parabolic`](@ref).
+
+The returned linear operator `A` is a matrix-free representation which can be
+supplied to iterative solvers from, e.g., [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
+
+It is also possible to use this to construct a sparse matrix without the detour of constructing
+first the full Jacobian by calling
+```julia
+using SparseArrays
+A_map, b = linear_structure_parabolic(semi, t0 = t0)
+A_sparse = sparse(A_map)
+```
+which can then be further used to construct for instance a
+[`MatrixOperator`](https://docs.sciml.ai/SciMLOperators/stable/tutorials/getting_started/#Simplest-Operator:-MatrixOperator)
+from [SciMLOperators.jl](https://docs.sciml.ai/SciMLOperators/stable/).
+This is especially useful for IMEX schemes where the parabolic part is implicitly,
+and for a `MatrixOperator` only factorized once.
+"""
+function linear_structure_parabolic(semi::SemidiscretizationHyperbolicParabolic;
+                                    t0 = zero(real(semi)))
+    if have_constant_diffusivity(semi.equations_parabolic) == False()
+        throw(ArgumentError("`linear_structure_parabolic` expects equations with constant diffusive terms."))
+    end
+
+    # allocate memory
+    u_ode = allocate_coefficients(mesh_equations_solver_cache(semi)...)
+    du_ode = similar(u_ode)
+
+    # get the parabolic right hand side from boundary conditions and optional source terms
+    u_ode .= zero(eltype(u_ode))
+    rhs_parabolic!(du_ode, u_ode, semi, t0)
+    b = -du_ode
+
+    # Create a copy of `b` used internally to extract the linear part of `semi`.
+    # This is necessary to get everything correct when the user updates the
+    # returned vector `b`.
+    b_tmp = copy(b)
+
+    # wrap the linear operator
+    A = LinearMap(length(u_ode), ismutating = true) do dest, src
+        rhs_parabolic!(dest, src, semi, t0)
+
+        @. dest += b_tmp
+        return dest
+    end
+
+    return A, b
+end
+
 """
     jacobian_ad_forward_parabolic(semi::SemidiscretizationHyperbolicParabolic;
                                   t0 = zero(real(semi)),

test/Project.toml

@@ -30,6 +30,8 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
@@ -38,40 +40,42 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26"
 
 [compat]
-ADTypes = "1.11"
-Adapt = "4"
-Aqua = "0.8"
-CUDA = "5.8"
-CairoMakie = "0.12, 0.13, 0.14, 0.15"
+ADTypes = "1.20"
+Adapt = "4.4"
+Aqua = "0.8.14"
+CUDA = "5.9"
+CairoMakie = "0.15.7"
 Convex = "0.16"
-DelimitedFiles = "1"
-DoubleFloats = "1.4.0"
-Downloads = "1"
+DelimitedFiles = "1.9.1"
+DoubleFloats = "1.5"
+Downloads = "1.6"
 ECOS = "1.1.2"
-ExplicitImports = "1.0.1"
-FiniteDiff = "2.27.0"
-ForwardDiff = "0.10.36, 1"
-Krylov = "0.9.10, 0.10"
+ExplicitImports = "1.14"
+FiniteDiff = "2.29"
+ForwardDiff = "1.3"
+Krylov = "0.10.3"
 LinearAlgebra = "1"
-LinearSolve = "2.36.1, 3"
-MPI = "0.20.6"
+LinearSolve = "3.48"
+MPI = "0.20.23"
 NLsolve = "4.5.1"
-OrdinaryDiffEqBDF = "1.1"
-OrdinaryDiffEqCore = "1.6"
-OrdinaryDiffEqFeagin = "1"
-OrdinaryDiffEqHighOrderRK = "1.1"
-OrdinaryDiffEqLowOrderRK = "1.2"
-OrdinaryDiffEqLowStorageRK = "1.2"
-OrdinaryDiffEqSDIRK = "1.1"
-OrdinaryDiffEqSSPRK = "1.2"
-OrdinaryDiffEqTsit5 = "1.1"
-Plots = "1.26"
+OrdinaryDiffEqBDF = "1.10.1"
+OrdinaryDiffEqCore = "1.36"
+OrdinaryDiffEqFeagin = "1.4"
+OrdinaryDiffEqHighOrderRK = "1.5"
+OrdinaryDiffEqLowOrderRK = "1.6"
+OrdinaryDiffEqLowStorageRK = "1.7"
+OrdinaryDiffEqSDIRK = "1.7"
+OrdinaryDiffEqSSPRK = "1.7"
+OrdinaryDiffEqTsit5 = "1.5"
+Plots = "1.41.2"
 Printf = "1"
-Quadmath = "0.5.10"
+Quadmath = "0.5.13"
 Random = "1"
-SparseArrays = "1"
+SciMLBase = "2.127.0"
+SciMLOperators = "1.13.0"
+SparseArrays = "1.10"
 SparseConnectivityTracer = "1.0.1"
-SparseMatrixColorings = "0.4.21"
-StableRNGs = "1.0.2"
+SparseMatrixColorings = "0.4.23"
+StableRNGs = "1.0.4"
 Test = "1"
 TrixiTest = "0.2.1"

test/test_parabolic_2d.jl

@@ -253,6 +253,15 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "TreeMesh2D: elixir_advection_diffusion_imex_operator.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
+                                 "elixir_advection_diffusion_imex_operator.jl"),
+                        l2=[7.542670562162156e-8], linf=[3.972014046560446e-7])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "TreeMesh2D: elixir_advection_diffusion_nonperiodic.jl (LDG)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem",
                                  "elixir_advection_diffusion_nonperiodic.jl"),