diff --git a/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl b/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl new file mode 100644 index 00000000000..c5998ae8d6e --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_euler_convergence_implicit_sparse_jacobian.jl @@ -0,0 +1,166 @@ +using Trixi +using OrdinaryDiffEqSDIRK + +# Functionality for automatic sparsity detection +using SparseDiffTools, Symbolics + +############################################################################################### +### Overloads to construct the `LobattoLegendreBasis` with `Real` type (supertype of `Num`) ### + +# Required for setting up the Lobatto-Legendre basis for abstract `Real` type. +# Constructing the Lobatto-Legendre basis with `Real` instead of `Num` is +# significantly easier as we do not have to care about e.g. if-clauses. +# As a consequence, we need to provide some overloads hinting towards the intended behavior. + +const float_type = Float64 # Actual floating point type for the simulation + +# Newton tolerance for finding LGL nodes & weights +Trixi.eps(::Type{Real}) = Base.eps(float_type) +# There are some places where `one(RealT)` or `zero(uEltype)` is called where `RealT` or `uEltype` is `Real`. +# This returns an `Int64`, i.e., `1` or `0`, respectively which gives errors when a floating-point alike type is expected. +Trixi.one(::Type{Real}) = Base.one(float_type) +Trixi.zero(::Type{Real}) = Base.zero(float_type) + +module RealMatMulOverload + +# Multiplying two Matrix{Real}s gives a Matrix{Any}. +# This causes problems when instantiating the Legendre basis, which calls +# `calc_{forward,reverse}_{upper, lower}` which in turn uses the matrix multiplication +# which is overloaded here in construction of the interpolation/projection operators +# required for mortars. +function Base.:*(A::Matrix{Real}, B::Matrix{Real})::Matrix{Real} + m, n = size(A, 1), size(B, 2) + kA = size(A, 2) + kB = size(B, 1) + @assert kA==kB "Matrix dimensions must match for multiplication" + + C = Matrix{Real}(undef, m, n) + for i in 1:m, j in 1:n + #acc::Real = zero(promote_type(typeof(A[i,1]), typeof(B[1,j]))) + acc = zero(Real) + for k in 1:kA + acc += A[i, k] * B[k, j] + end + C[i, j] = acc + end + return C +end +end + +import .RealMatMulOverload + +# We need to avoid if-clauses to be able to use `Num` type from Symbolics without additional hassle. +# In the Trixi implementation, we overload the sqrt function to first check if the argument +# is < 0 and then return NaN instead of an error. +# To turn off this behaviour, we switch back to the Base implementation here which does not contain an if-clause. +Trixi.sqrt(x::Num) = Base.sqrt(x) + +############################################################################################### +### equations and solver ### + +equations = CompressibleEulerEquations2D(1.4) + +# For sparsity detection we can only use `flux_lax_friedrichs` at the moment since this is +# `if`-clause free +surface_flux = flux_lax_friedrichs + +# `RealT = Real` requires fewer overloads than the more explicit `RealT = Num` from Symbolics. +# `solver_real` is used for computing the Jacobian sparsity pattern +solver_real = DGSEM(polydeg = 3, surface_flux = surface_flux, RealT = Real) +# `solver_float` is used for the subsequent simulation +solver_float = DGSEM(polydeg = 3, surface_flux = surface_flux, RealT = float_type) + +############################################################################################### +### mesh ### + +# Mapping as described in https://arxiv.org/abs/2012.12040, +# reduced to 2D on [0, 2]^2 instead of [0, 3]^3 +function mapping(xi_, eta_) + # Transform input variables between -1 and 1 onto [0,2] + xi = xi_ + 1 + eta = eta_ + 1 + + y = eta + 1 / 4 * (cos(pi * (xi - 1)) * + cos(0.5 * pi * (eta - 1))) + + x = xi + 1 / 4 * (cos(0.5 * pi * (xi - 1)) * + cos(2 * pi * (y - 1))) + + return SVector(x, y) +end +cells_per_dimension = (16, 16) +mesh = StructuredMesh(cells_per_dimension, mapping) + +############################################################################################### +### semidiscretizations ### + +initial_condition = initial_condition_convergence_test + +# `semi_real` is used for computing the Jacobian sparsity pattern +semi_real = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_convergence_test, + solver_real, + source_terms = source_terms_convergence_test) +# `semi_float` is used for the subsequent simulation +semi_float = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_convergence_test, + solver_float, + source_terms = source_terms_convergence_test) + +t0 = 0.0 # Re-used for the ODE function defined below +t_end = 5.0 +t_span = (t0, t_end) + +# Call `semidiscretize` on `semi_float` to create the ODE problem to have access to the initial condition. +ode_float = semidiscretize(semi_float, t_span) +u0_ode = ode_float.u0 +du_ode = similar(u0_ode) + +############################################################################################### +### Compute the Jacobian with SparseDiffTools ### + +# Create a function with two parameters: `du_ode` and `u0_ode` +# to fulfill the requirements of an in_place function in SparseDiffTools +# (see example function `f` from https://docs.sciml.ai/SparseDiffTools/dev/#Example) +rhs = (du_ode, u0_ode) -> Trixi.rhs!(du_ode, u0_ode, semi_real, t0) + +# Taken from example linked above to detect the pattern and choose how to do the differentiation +sd = SymbolicsSparsityDetection() +ad_type = AutoForwardDiff() +sparse_adtype = AutoSparse(ad_type) + +# `sparse_cache` will reduce calculation time when Jacobian is calculated multiple times +sparse_cache = sparse_jacobian_cache(sparse_adtype, sd, rhs, du_ode, u0_ode) +# The jacobian can eventually be re-computed in-place using +#J_prealloc = sparse_jacobian(sparse_adtype, sparse_cache, rhs, du_ode, u0_ode) +#sparse_jacobian!(J_prealloc, sparse_adtype, sparse_cache, rhs, du_ode, u0_ode) + +############################################################################################### +### Set up sparse-aware ODEProblem ### + +# Revert overrides from above for the actual simulation - +# not strictly necessary, but good practice +Trixi.eps(x::Type{Real}) = Base.eps(x) +Trixi.one(x::Type{Real}) = Base.one(x) +Trixi.zero(x::Type{Real}) = Base.zero(x) + +# Supply Jacobian prototype and coloring vector to the semidiscretization +ode_float_jac_sparse = semidiscretize(semi_float, t_span, + jac_prototype = sparse_cache.jac_prototype, + colorvec = sparse_cache.coloring.colorvec) + +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi_float, interval = 50) +alive_callback = AliveCallback(alive_interval = 3) + +# Note: No `stepsize_callback` due to implicit solver +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback) + +############################################################################################### +### solve the ODE problem ### + +sol = solve(ode_float_jac_sparse, # using `ode_float` is essentially infeasible, even single step takes ages! + # Default `AutoForwardDiff()` is not yet working, + # probably related to https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers + Kvaerno4(; autodiff = AutoFiniteDiff()); + dt = 0.05, save_everystep = false, callback = callbacks); diff --git a/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl new file mode 100644 index 00000000000..8e17a94dabf --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian.jl @@ -0,0 +1,147 @@ +using Trixi +using OrdinaryDiffEqSDIRK + +# Functionality for automatic sparsity detection +using SparseDiffTools, Symbolics + +############################################################################################### +### Overloads to construct the `LobattoLegendreBasis` with `Real` type (supertype of `Num`) ### + +# Required for setting up the Lobatto-Legendre basis for abstract `Real` type. +# Constructing the Lobatto-Legendre basis with `Real` instead of `Num` is +# significantly easier as we do not have to care about e.g. if-clauses. +# As a consequence, we need to provide some overloads hinting towards the intended behavior. + +const float_type = Float64 # Actual floating point type for the simulation + +# Newton tolerance for finding LGL nodes & weights +Trixi.eps(::Type{Real}) = Base.eps(float_type) +# There are some places where `one(RealT)` or `zero(uEltype)` is called where `RealT` or `uEltype` is `Real`. +# This returns an `Int64`, i.e., `1` or `0`, respectively which gives errors when a floating-point alike type is expected. +Trixi.one(::Type{Real}) = Base.one(float_type) +Trixi.zero(::Type{Real}) = Base.zero(float_type) + +module RealMatMulOverload + +# Multiplying two Matrix{Real}s gives a Matrix{Any}. +# This causes problems when instantiating the Legendre basis, which calls +# `calc_{forward,reverse}_{upper, lower}` which in turn uses the matrix multiplication +# which is overloaded here in construction of the interpolation/projection operators +# required for mortars. +function Base.:*(A::Matrix{Real}, B::Matrix{Real})::Matrix{Real} + m, n = size(A, 1), size(B, 2) + kA = size(A, 2) + kB = size(B, 1) + @assert kA==kB "Matrix dimensions must match for multiplication" + + C = Matrix{Real}(undef, m, n) + for i in 1:m, j in 1:n + #acc::Real = zero(promote_type(typeof(A[i,1]), typeof(B[1,j]))) + acc = zero(Real) + for k in 1:kA + acc += A[i, k] * B[k, j] + end + C[i, j] = acc + end + return C +end +end + +import .RealMatMulOverload + +############################################################################################### +### semidiscretizations of the linear advection equation ### + +advection_velocity = (0.2, -0.7) +equation = LinearScalarAdvectionEquation2D(advection_velocity) + +# `RealT = Real` requires fewer overloads than the more explicit `RealT = Num` from Symbolics +# `solver_real` is used for computing the Jacobian sparsity pattern +solver_real = DGSEM(polydeg = 3, surface_flux = flux_godunov, RealT = Real) +# `solver_float` is used for the subsequent simulation +solver_float = DGSEM(polydeg = 3, surface_flux = flux_godunov) + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) + +# `semi_real` is used for computing the Jacobian sparsity pattern +semi_real = SemidiscretizationHyperbolic(mesh, equation, + initial_condition_convergence_test, + solver_real) +# `semi_float` is used for the subsequent simulation +semi_float = SemidiscretizationHyperbolic(mesh, equation, + initial_condition_convergence_test, + solver_float) + +t0 = 0.0 # Re-used for the ODE function defined below +t_end = 1.0 +t_span = (t0, t_end) + +# Call `semidiscretize` on `semi_float` to create the ODE problem to have access to the initial condition. +# For the linear example considered here one could also use an arbitrary vector for the initial condition. +ode_float = semidiscretize(semi_float, t_span) +u0_ode = ode_float.u0 +du_ode = similar(u0_ode) + +############################################################################################### +### Compute the Jacobian with SparseDiffTools ### + +# Create a function with two parameters: `du_ode` and `u0_ode` +# to fulfill the requirements of an in_place function in SparseDiffTools +# (see example function `f` from https://docs.sciml.ai/SparseDiffTools/dev/#Example) +rhs = (du_ode, u0_ode) -> Trixi.rhs!(du_ode, u0_ode, semi_real, t0) + +# Taken from example linked above to detect the pattern and choose how to do the AutoDiff automatically +sd = SymbolicsSparsityDetection() +ad_type = AutoFiniteDiff() # `AutoForwardDiff()` does also work, but not strictly necessary for sparsity detection only +sparse_adtype = AutoSparse(ad_type) + +# `sparse_cache` will reduce calculation time when Jacobian is calculated multiple times, +# which is in principle not required for the linear problem considered here. +sparse_cache = sparse_jacobian_cache(sparse_adtype, sd, rhs, du_ode, u0_ode) + +############################################################################################### +### Set up sparse-aware ODEProblem ### + +# Revert overrides from above for the actual simulation - +# not strictly necessary, but good practice +Trixi.eps(x::Type{Real}) = Base.eps(x) +Trixi.one(x::Type{Real}) = Base.one(x) +Trixi.zero(x::Type{Real}) = Base.zero(x) + +# Supply Jacobian prototype and coloring vector to the semidiscretization +ode_float_jac_sparse = semidiscretize(semi_float, t_span, + jac_prototype = sparse_cache.jac_prototype, + colorvec = sparse_cache.coloring.colorvec) + +# Note: We tried for this linear problem providing the constant, sparse Jacobian directly via +# +# const jac_sparse = sparse_jacobian(sparse_adtype, sparse_cache, rhs, du_ode, u0_ode) +# function jac_sparse_func!(J, u, p, t) +# J = jac_sparse +# end +# SciMLBase.ODEFunction(rhs!, jac_prototype=float.(jac_prototype), colorvec=colorvec, jac = jac_sparse_func!) +# +# which turned out (for the considered problems) to be significantly slower than +# just using the prototype and the coloring vector. + +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi_float, interval = 10) +save_restart = SaveRestartCallback(interval = 100, + save_final_restart = true) + +# Note: No `stepsize_callback` due to (implicit) solver with adaptive timestep control +callbacks = CallbackSet(summary_callback, analysis_callback, save_restart) + +############################################################################################### +### solve the ODE problem ### + +sol = solve(ode_float_jac_sparse, # using `ode_float_jac_sparse` instead of `ode_float` results in speedup of factors 10-15! + # Default `AutoForwardDiff()` is not yet working, + # probably related to https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers + TRBDF2(; autodiff = ad_type); + dt = 0.1, save_everystep = false, callback = callbacks); diff --git a/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian_restart.jl b/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian_restart.jl new file mode 100644 index 00000000000..6c61c80a6d7 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_advection_implicit_sparse_jacobian_restart.jl @@ -0,0 +1,28 @@ +using Trixi + +############################################################################### +# create a restart file + +elixir_file = "elixir_advection_implicit_sparse_jacobian.jl" +restart_file = "restart_000000006.h5" + +trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file)) + +############################################################################### + +restart_filename = joinpath("out", restart_file) +t_span = (load_time(restart_filename), 2.0) +dt_restart = load_dt(restart_filename) + +ode_float_jac_sparse = semidiscretize(semi_float, t_span, + jac_prototype = sparse_cache.jac_prototype, + colorvec = sparse_cache.coloring.colorvec) + +############################################################################### +# run the simulation + +sol = solve(ode_float_jac_sparse, # using `ode_float_jac_sparse` instead of `ode_float` results in speedup of factors 10-15! + # Default `AutoForwardDiff()` is not yet working, + # probably related to https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers + TRBDF2(; autodiff = ad_type); + adaptive = true, dt = dt_restart, save_everystep = false, callback = callbacks); diff --git a/src/semidiscretization/semidiscretization.jl b/src/semidiscretization/semidiscretization.jl index 34e8e609ada..73d4cc657ab 100644 --- a/src/semidiscretization/semidiscretization.jl +++ b/src/semidiscretization/semidiscretization.jl @@ -78,16 +78,25 @@ function calc_error_norms(u_ode, t, analyzer, semi::AbstractSemidiscretization, end """ - semidiscretize(semi::AbstractSemidiscretization, tspan) + semidiscretize(semi::AbstractSemidiscretization, tspan; + jac_prototype::Union{AbstractMatrix, Nothing} = nothing, + colorvec::Union{AbstractVector, Nothing} = nothing, + storage_type = nothing, + real_type = nothing) Wrap the semidiscretization `semi` as an ODE problem in the time interval `tspan` that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). -The optional keyword arguments `storage_type` and `real_type` configure the underlying computational -datastructures. `storage_type` changes the fundamental array type being used, allowing the -experimental use of `CuArray` or other GPU array types. `real_type` changes the computational data type being used. +Optional keyword arguments: +- `jac_prototype` and `colorvec`: Expected to come from [SparseDiffTools.jl](https://github.com/JuliaDiff/SparseDiffTools.jl) + and specify the sparsity structure of the Jacobian to enable efficient implicit time stepping. +- `storage_type` and `real_type`: Configure the underlying computational datastructures. + `storage_type` changes the fundamental array type being used, allowing the experimental use of `CuArray` + or other GPU array types. `real_type` changes the computational data type being used. """ function semidiscretize(semi::AbstractSemidiscretization, tspan; + jac_prototype::Union{AbstractMatrix, Nothing} = nothing, + colorvec::Union{AbstractVector, Nothing} = nothing, reset_threads = true, storage_type = nothing, real_type = nothing) @@ -111,25 +120,47 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan; end end - u0_ode = compute_coefficients(first(tspan), semi) + u0_ode = compute_coefficients(first(tspan), semi) # Invoke initial condition + # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using # mpi_isparallel() && MPI.Barrier(mpi_comm()) # See https://github.com/trixi-framework/Trixi.jl/issues/328 iip = true # is-inplace, i.e., we modify a vector when calling rhs! specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi) - return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) + + # Check if Jacobian prototype and coloring vector are provided for sparse Jacobian + if jac_prototype !== nothing && colorvec !== nothing + # See SparseDiffTools.jl docs (https://github.com/JuliaDiff/SparseDiffTools.jl) for documentation of `jac_prototype` and `colorvec` + + # Convert the `jac_prototype` to real type, as seen here: + # https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Declaring-a-Sparse-Jacobian-with-Automatic-Sparsity-Detection + ode = SciMLBase.ODEFunction(rhs!, jac_prototype = float.(jac_prototype), + colorvec = colorvec) + + return ODEProblem{iip, specialize}(ode, u0_ode, tspan, semi) + else + return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) + end end """ semidiscretize(semi::AbstractSemidiscretization, tspan, - restart_file::AbstractString) + restart_file::AbstractString; + jac_prototype::Union{AbstractMatrix, Nothing} = nothing, + colorvec::Union{AbstractVector, Nothing} = nothing) Wrap the semidiscretization `semi` as an ODE problem in the time interval `tspan` that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). The initial condition etc. is taken from the `restart_file`. + +Optional keyword arguments: +- `jac_prototype` and `colorvec`: Expected to come from [SparseDiffTools.jl](https://github.com/JuliaDiff/SparseDiffTools.jl) + and specify the sparsity structure of the Jacobian to enable efficient implicit time stepping. """ function semidiscretize(semi::AbstractSemidiscretization, tspan, restart_file::AbstractString; + jac_prototype::Union{AbstractMatrix, Nothing} = nothing, + colorvec::Union{AbstractVector, Nothing} = nothing, reset_threads = true) # Optionally reset Polyester.jl threads. See # https://github.com/trixi-framework/Trixi.jl/issues/1583 @@ -138,13 +169,27 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan, Polyester.reset_threads!() end - u0_ode = load_restart_file(semi, restart_file) + u0_ode = load_restart_file(semi, restart_file) # Load initial condition from restart file + # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using # mpi_isparallel() && MPI.Barrier(mpi_comm()) # See https://github.com/trixi-framework/Trixi.jl/issues/328 iip = true # is-inplace, i.e., we modify a vector when calling rhs! specialize = SciMLBase.FullSpecialize # specialize on rhs! and parameters (semi) - return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) + + # Check if Jacobian prototype and coloring vector are provided for sparse Jacobian + if jac_prototype !== nothing && colorvec !== nothing + # See SparseDiffTools.jl docs (https://github.com/JuliaDiff/SparseDiffTools.jl) for documentation of `jac_prototype` and `colorvec` + + # Convert the `jac_prototype` to real type, as seen here: + # https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Declaring-a-Sparse-Jacobian-with-Automatic-Sparsity-Detection + ode = SciMLBase.ODEFunction(rhs!, jac_prototype = float.(jac_prototype), + colorvec = colorvec) + + return ODEProblem{iip, specialize}(ode, u0_ode, tspan, semi) + else + return ODEProblem{iip, specialize}(rhs!, u0_ode, tspan, semi) + end end """ diff --git a/test/Project.toml b/test/Project.toml index c996ec756a5..572845d28e9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -25,7 +25,10 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26" @@ -56,6 +59,9 @@ Plots = "1.26" Printf = "1" Quadmath = "0.5.10" Random = "1" +SparseArrays = "1" +SparseDiffTools = "2.26.0" StableRNGs = "1.0.2" +Symbolics = "6" Test = "1" TrixiTest = "0.1" diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 00b9825812f..d824872427a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -319,6 +319,32 @@ end end end +@trixi_testset "elixir_euler_convergence_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_convergence_implicit_sparse_jacobian.jl"), + t_end=1.0, + l2=[ + 0.0025545032994391684, + 0.002584889213509444, + 0.002585815262287157, + 0.0031668773337868226 + ], + linf=[ + 0.010367159504632184, + 0.00932621263313771, + 0.008372785091579793, + 0.011242647117395421 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi_float, t)) < 1000 + end +end + @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"), l2=[ diff --git a/test/test_tree_2d_advection.jl b/test/test_tree_2d_advection.jl index 78b20e0036d..5cf75e23e23 100644 --- a/test/test_tree_2d_advection.jl +++ b/test/test_tree_2d_advection.jl @@ -40,6 +40,34 @@ end end end +@trixi_testset "elixir_advection_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_implicit_sparse_jacobian.jl"), + l2=[0.003003253325111022], linf=[0.004256250998163846]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi_float, t)) < 1000 + end +end + +@trixi_testset "elixir_advection_implicit_sparse_jacobian_restart.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_implicit_sparse_jacobian_restart.jl"), + l2=[0.004964299251933375], linf=[0.007028151045231468]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi_float, t)) < 1000 + end +end + @trixi_testset "elixir_advection_restart.jl" begin using OrdinaryDiffEqSSPRK: SSPRK43 println("═"^100) diff --git a/test/test_trixi.jl b/test/test_trixi.jl index 00baf2162d0..6e359677186 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -121,8 +121,8 @@ macro test_nowarn_mod(expr, additional_ignore_content = []) r"┌ Warning: Problem status ALMOST_INFEASIBLE; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n", r"┌ Warning: Problem status ALMOST_OPTIMAL; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n", # Warnings for higher-precision floating data types - r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/interpolation.jl:118 =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/.*\n", - r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/interpolation.jl:136 =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/.*\n" + r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/.* =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/.*\n", + r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/.* =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/.*\n" ] append!($additional_ignore_content, add_to_additional_ignore_content) @trixi_test_nowarn $(esc(expr)) $additional_ignore_content diff --git a/test/test_unit.jl b/test/test_unit.jl index 0a305a3f0c4..b1440d34c34 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -4,6 +4,7 @@ using Test using Trixi using LinearAlgebra: norm, dot +using SparseArrays using DelimitedFiles: readdlm # Use Convex and ECOS to load the extension that extends functions for testing @@ -15,6 +16,9 @@ using ECOS: Optimizer # PERK Single p3 Constructors using NLsolve: nlsolve +# For sparsity detection with Symbolics +using SparseDiffTools, Symbolics + include("test_trixi.jl") # Start with a clean environment: remove Trixi.jl output directory if it exists @@ -2604,5 +2608,113 @@ end equations) end end + +# Multiplying two Matrix{Real}s gives a Matrix{Any}. +# This causes problems when instantiating the Legendre basis, which calls +# `calc_{forward,reverse}_{upper, lower}` which in turn uses the matrix multiplication +# which is overloaded here in construction of the interpolation/projection operators +# required for mortars. +function Base.:*(A::Matrix{Real}, B::Matrix{Real})::Matrix{Real} + m, n = size(A, 1), size(B, 2) + kA = size(A, 2) + kB = size(B, 1) + @assert kA==kB "Matrix dimensions must match for multiplication" + + C = Matrix{Real}(undef, m, n) + for i in 1:m, j in 1:n + #acc::Real = zero(promote_type(typeof(A[i,1]), typeof(B[1,j]))) + acc = zero(Real) + for k in 1:kA + acc += A[i, k] * B[k, j] + end + C[i, j] = acc + end + return C +end + +@testset "SparseDiff Jacobian = {ForwardDiff Jacobian, LinearStructure}" begin + ############################################################################################### + ### Overloads to construct the `LobattoLegendreBasis` with `Real` type (supertype of `Num`) ### + + # Required for setting up the Lobatto-Legendre basis for abstract `Real` type. + # Constructing the Lobatto-Legendre basis with `Real` instead of `Num` is + # significantly easier as we do not have to care about e.g. if-clauses. + # As a consequence, we need to provide some overloads hinting towards the intended behavior. + + float_type = Float64 # Actual floating point type for the simulation + + # Newton tolerance for finding LGL nodes & weights + Trixi.eps(::Type{Real}) = Base.eps(float_type) + # There are some places where `one(RealT)` or `zero(uEltype)` is called where `RealT` or `uEltype` is `Real`. + # This returns an `Int64`, i.e., `1` or `0`, respectively which gives errors when a floating-point alike type is expected. + Trixi.one(::Type{Real}) = Base.one(float_type) + Trixi.zero(::Type{Real}) = Base.zero(float_type) + + ############################################################################################### + + advection_velocities = (0.2, -0.7) + equations = LinearScalarAdvectionEquation2D(advection_velocities) + + # `RealT = Real` requires fewer overloads than the more explicit `RealT = Num` from Symbolics. + # `solver_real` is used for computing the Jacobian sparsity pattern + solver_real = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + RealT = Real) + # `solver_float` is used for the subsequent simulation + solver_float = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + RealT = float_type) + + coordinates_min = (-1.0, -1.0) + coordinates_max = (1.0, 1.0) + + mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) + + # `semi_real` is used for computing the Jacobian sparsity pattern + semi_real = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_convergence_test, + solver_real) + # `semi_float` is used for the subsequent simulation + semi_float = SemidiscretizationHyperbolic(mesh, equations, + initial_condition_convergence_test, + solver_float) + + t0 = 0.0 # Re-used for the ODE function defined below + t_end = 5.0 + t_span = (t0, t_end) + + # Call `semidiscretize` on `semi_float` to create the ODE problem to have access to the initial condition. + ode_float = semidiscretize(semi_float, t_span) + u0_ode = ode_float.u0 + du_ode = similar(u0_ode) + + ############################################################################################### + ### Compute the Jacobian with SparseDiffTools ### + + # Create a function with two parameters: `du_ode` and `u0_ode` + # to fulfill the requirements of an in_place function in SparseDiffTools + # (see example function `f` from https://docs.sciml.ai/SparseDiffTools/dev/#Example) + rhs = (du_ode, u0_ode) -> Trixi.rhs!(du_ode, u0_ode, semi_real, t0) + + # Taken from example linked above to detect the pattern and choose how to do the differentiation + ad_type = AutoForwardDiff() + sparse_adtype = AutoSparse(ad_type) + sd = SymbolicsSparsityDetection() + sparse_cache = sparse_jacobian_cache(sparse_adtype, sd, rhs, du_ode, u0_ode) + jac_sparse = sparse_jacobian(sparse_adtype, sparse_cache, rhs, du_ode, u0_ode) + + jac_forward_diff = jacobian_ad_forward(semi_float) + + @test jac_sparse == jac_forward_diff + @test Matrix(jac_sparse) == jac_forward_diff + @test jac_sparse == sparse(jac_forward_diff) + + A, _ = linear_structure(semi_float) + + @test jac_sparse == Matrix(A) + @test Matrix(jac_sparse) == Matrix(A) + @test jac_sparse == sparse(A) +end end + end #module