From 7097ac279251adde9da2c5654787c0fb070757ea Mon Sep 17 00:00:00 2001 From: Steve Date: Wed, 3 Sep 2025 21:27:36 +0200 Subject: [PATCH 01/22] Hyperbolic-Parabolic elixir using SCT --- ...tion_diffusion_implicit_sparse_jacobian.jl | 138 ++++++++++++++++++ ...semidiscretization_hyperbolic_parabolic.jl | 78 ++++++++-- 2 files changed, 206 insertions(+), 10 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl new file mode 100644 index 00000000000..9e061603631 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -0,0 +1,138 @@ +using Trixi +using SparseConnectivityTracer # For obtaining the Jacobian sparsity pattern +using SparseMatrixColorings # For obtaining the coloring vector +using OrdinaryDiffEqSDIRK, ADTypes + +############################################################################### +# semidiscretization of the linear advection-diffusion equation + +advection_velocity = (1.5, 1.0) +equations = LinearScalarAdvectionEquation2D(advection_velocity) +diffusivity() = 5.0e-2 +equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +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)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + periodicity = true, + n_cells_max = 30_000) # set maximum capacity of tree data structure + +# Define initial condition +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 + +# define periodic boundary conditions everywhere +boundary_conditions = boundary_condition_periodic +boundary_conditions_parabolic = boundary_condition_periodic + +############################################################################### +### semidiscretization for sparsity detection ### + +jac_detector = TracerSparsityDetector() +# We need to construct the semidiscretization with the correct +# sparsity-detection ready datatype, which is retrieved here +jac_eltype = jacobian_eltype(real(solver), jac_detector) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver, + uEltype = jac_eltype; + solver_parabolic = ViscousFormulationBassiRebay1(), + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +tspan = (0.0, 1.5) # Re-used for wrapping `rhs_parabolic!` below + +# Call `semidiscretize` to create the ODE problem to have access to the +# initial condition based on which the sparsity pattern is computed +ode_jac_type = semidiscretize(semi_jac_type, tspan) +u0_ode = ode_jac_type.u0 +du_ode = similar(u0_ode) + +############################################################################### +### Compute the Jacobian sparsity pattern ### + +# Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian +# https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitFunction + +# Wrap the `Trixi.rhs_parabolic!` function to match the signature `f!(du, u)`, see +# https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity +rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type, tspan[1]) + +jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) + +# For most efficient solving we also want the coloring vector + +coloring_prob = ColoringProblem(; structure = :nonsymmetric, partition = :column) +coloring_alg = GreedyColoringAlgorithm(; decompression = :direct) +coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) +coloring_vec_parabolic = column_colors(coloring_result) + +############################################################################### +### sparsity-aware semidiscretization and ode ### + +# Semidiscretization for actual simulation. `eEltype` is here retrieved from `solver` +semi_float_type = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver; + solver_parabolic = ViscousFormulationBassiRebay1(), + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +# Supply Jacobian prototype and coloring vector to the semidiscretization +ode_jac_sparse = semidiscretize(semi_float_type, tspan, + jac_prototype_parabolic = jac_prototype_parabolic, + colorvec_parabolic = coloring_vec_parabolic) +# using "dense" `ode = semidiscretize(semi_float_type, tspan)` is 4-6 times slower! + +############################################################################### +### callbacks ### + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi_float_type, interval = analysis_interval) + +# The AliveCallback prints short status information in regular intervals +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +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, alive_callback, save_restart) + +############################################################################### +### solve the ODE problem ### + +time_int_tol = 1.0e-9 +time_abs_tol = 1.0e-9 + +sol = solve(ode_jac_sparse, TRBDF2(; autodiff = AutoFiniteDiff()); + dt = 0.1, save_everystep = false, + abstol = time_abs_tol, reltol = time_int_tol, + ode_default_options()..., callback = callbacks) \ No newline at end of file diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 54ede387fa2..91e55f9675e 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -264,15 +264,25 @@ function get_node_variables!(node_variables, u_ode, end """ - semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan) + semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan; + jac_prototype_parabolic::Union{AbstractMatrix, Nothing} = nothing, + colorvec_parabolic::Union{AbstractVector, Nothing} = nothing) Wrap the semidiscretization `semi` as a split ODE problem in the time interval `tspan` that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). The parabolic right-hand side is the first function of the split ODE problem and will be used by default by the implicit part of IMEX methods from the SciML ecosystem. + +Optional keyword arguments: +- `jac_prototype_parabolic`: Expected to come from [SparseConnectivityTracer.jl](https://github.com/adrhill/SparseConnectivityTracer.jl). + Specifies the sparsity structure of the parabolic function's Jacobian to enable e.g. efficient implicit time stepping. +- `colorvec_parabolic`: Expected to come from [SparseMatrixColorings.jl](https://github.com/gdalle/SparseMatrixColorings.jl). + Allows for even faster Jacobian computation. Not necessarily required when `jac_prototype_parabolic` is given. """ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan; + jac_prototype_parabolic::Union{AbstractMatrix, Nothing} = nothing, + colorvec_parabolic::Union{AbstractVector, Nothing} = nothing, reset_threads = true) # Optionally reset Polyester.jl threads. See # https://github.com/trixi-framework/Trixi.jl/issues/1583 @@ -286,15 +296,36 @@ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan; # 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_parabolic!, rhs! - # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the - # first function implicitly and the second one explicitly. Thus, we pass the - # stiffer parabolic function first. - return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) + + # Check if Jacobian prototype is provided for sparse Jacobian + if jac_prototype_parabolic !== nothing + # Convert `jac_prototype_parabolic` to real type, as seen here: + # https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Declaring-a-Sparse-Jacobian-with-Automatic-Sparsity-Detection + parabolic_ode = SciMLBase.ODEFunction(rhs_parabolic!, + jac_prototype = convert.(eltype(u0_ode), + jac_prototype_parabolic), + colorvec = colorvec_parabolic) # coloring vector is optional + + # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the + # first function implicitly and the second one explicitly. Thus, we pass the + # stiffer parabolic function first. + return SplitODEProblem{iip}(parabolic_ode, rhs!, u0_ode, tspan, semi) + else + # We could also construct an `ODEFunction` without the Jacobian here, + # but we stick to the more light-weight direct in-place function `rhs_parabolic!`. + + # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the + # first function implicitly and the second one explicitly. Thus, we pass the + # stiffer parabolic function first. + return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) + end end """ semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan, - restart_file::AbstractString) + restart_file::AbstractString; + jac_prototype_parabolic::Union{AbstractMatrix, Nothing} = nothing, + colorvec_parabolic::Union{AbstractVector, Nothing} = nothing) Wrap the semidiscretization `semi` as a split ODE problem in the time interval `tspan` that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/). @@ -303,9 +334,17 @@ will be used by default by the implicit part of IMEX methods from the SciML ecosystem. The initial condition etc. is taken from the `restart_file`. + +Optional keyword arguments: +- `jac_prototype_parabolic`: Expected to come from [SparseConnectivityTracer.jl](https://github.com/adrhill/SparseConnectivityTracer.jl). + Specifies the sparsity structure of the parabolic function's Jacobian to enable e.g. efficient implicit time stepping. +- `colorvec_parabolic`: Expected to come from [SparseMatrixColorings.jl](https://github.com/gdalle/SparseMatrixColorings.jl). + Allows for even faster Jacobian computation. Not necessarily required when `jac_prototype_parabolic` is given. """ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan, restart_file::AbstractString; + jac_prototype_parabolic::Union{AbstractMatrix, Nothing} = nothing, + colorvec_parabolic::Union{AbstractVector, Nothing} = nothing, reset_threads = true) # Optionally reset Polyester.jl threads. See # https://github.com/trixi-framework/Trixi.jl/issues/1583 @@ -319,10 +358,29 @@ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan, # 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_parabolic!, rhs! - # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the - # first function implicitly and the second one explicitly. Thus, we pass the - # stiffer parabolic function first. - return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) + + # Check if Jacobian prototype is provided for sparse Jacobian + if jac_prototype_parabolic !== nothing + # Convert `jac_prototype_parabolic` to real type, as seen here: + # https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Declaring-a-Sparse-Jacobian-with-Automatic-Sparsity-Detection + parabolic_ode = SciMLBase.ODEFunction(rhs_parabolic!, + jac_prototype = convert.(eltype(u0_ode), + jac_prototype_parabolic), + colorvec = colorvec_parabolic) # coloring vector is optional + + # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the + # first function implicitly and the second one explicitly. Thus, we pass the + # stiffer parabolic function first. + return SplitODEProblem{iip}(parabolic_ode, rhs!, u0_ode, tspan, semi) + else + # We could also construct an `ODEFunction` without the Jacobian here, + # but we stick to the more light-weight direct in-place function `rhs_parabolic!`. + + # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the + # first function implicitly and the second one explicitly. Thus, we pass the + # stiffer parabolic function first. + return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) + end end function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t) From bcb5214f3f86b2a679aa66f9127e452944e741d5 Mon Sep 17 00:00:00 2001 From: Steve Date: Thu, 4 Sep 2025 19:13:42 +0200 Subject: [PATCH 02/22] Add CI test --- test/test_parabolic_2d.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 0d23b43ef4b..22d71b47362 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -264,6 +264,21 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "TreeMesh2D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", + "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), + initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, + l2=[4.3996570276332366e-5], linf=[7.656808122324943e-5]) + # 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_type, t)) < 1000 + end +end + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), From 6e87bd86b350324dd8f21b0381bba772699d8697 Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 9 Sep 2025 16:22:05 -0500 Subject: [PATCH 03/22] PR corrections --- ...ection_diffusion_implicit_sparse_jacobian.jl | 17 ++++------------- .../semidiscretization_hyperbolic_parabolic.jl | 3 +++ 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 9e061603631..fd1c5c953f4 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -11,19 +11,16 @@ equations = LinearScalarAdvectionEquation2D(advection_velocity) diffusivity() = 5.0e-2 equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux 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)) +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) -# Create a uniformly refined mesh with periodic boundaries mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 4, periodicity = true, - n_cells_max = 30_000) # set maximum capacity of tree data structure + n_cells_max = 30_000) -# Define initial condition function initial_condition_diffusive_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) # Store translated coordinate for easy use of exact solution @@ -41,7 +38,6 @@ function initial_condition_diffusive_convergence_test(x, t, end initial_condition = initial_condition_diffusive_convergence_test -# define periodic boundary conditions everywhere boundary_conditions = boundary_condition_periodic boundary_conditions_parabolic = boundary_condition_periodic @@ -53,7 +49,6 @@ jac_detector = TracerSparsityDetector() # sparsity-detection ready datatype, which is retrieved here jac_eltype = jacobian_eltype(real(solver), jac_detector) -# A semidiscretization collects data structures and functions for the spatial discretization semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver, @@ -90,7 +85,7 @@ coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) coloring_vec_parabolic = column_colors(coloring_result) ############################################################################### -### sparsity-aware semidiscretization and ode ### +### sparsity-aware semidiscretization and ODE ### # Semidiscretization for actual simulation. `eEltype` is here retrieved from `solver` semi_float_type = SemidiscretizationHyperbolicParabolic(mesh, @@ -109,15 +104,11 @@ ode_jac_sparse = semidiscretize(semi_float_type, tspan, ############################################################################### ### callbacks ### -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers summary_callback = SummaryCallback() -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results analysis_interval = 100 analysis_callback = AnalysisCallback(semi_float_type, interval = analysis_interval) -# The AliveCallback prints short status information in regular intervals alive_callback = AliveCallback(analysis_interval = analysis_interval) save_restart = SaveRestartCallback(interval = 100, diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 91e55f9675e..1e9144f3973 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -338,6 +338,9 @@ The initial condition etc. is taken from the `restart_file`. Optional keyword arguments: - `jac_prototype_parabolic`: Expected to come from [SparseConnectivityTracer.jl](https://github.com/adrhill/SparseConnectivityTracer.jl). Specifies the sparsity structure of the parabolic function's Jacobian to enable e.g. efficient implicit time stepping. + The [SplitODEProblem](https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitODEProblem) only expects the Jacobian + to be defined on the first function it takes in, which is treated implicitly. This corresponds to the parabolic right-hand side in Trixi. + The hyperbolic right-hand side is treated explicitly, and therefore its Jacobian is irrelevant. - `colorvec_parabolic`: Expected to come from [SparseMatrixColorings.jl](https://github.com/gdalle/SparseMatrixColorings.jl). Allows for even faster Jacobian computation. Not necessarily required when `jac_prototype_parabolic` is given. """ From 31b0394813392026f3781cf2b035c74d9b697536 Mon Sep 17 00:00:00 2001 From: Steve Date: Thu, 11 Sep 2025 13:17:36 -0500 Subject: [PATCH 04/22] Adding missed comment --- .../semidiscretization_hyperbolic_parabolic.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 1e9144f3973..052957af6cc 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -277,6 +277,9 @@ SciML ecosystem. Optional keyword arguments: - `jac_prototype_parabolic`: Expected to come from [SparseConnectivityTracer.jl](https://github.com/adrhill/SparseConnectivityTracer.jl). Specifies the sparsity structure of the parabolic function's Jacobian to enable e.g. efficient implicit time stepping. + The [SplitODEProblem](https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitODEProblem) only expects the Jacobian + to be defined on the first function it takes in, which is treated implicitly. This corresponds to the parabolic right-hand side in Trixi. + The hyperbolic right-hand side is treated explicitly, and therefore its Jacobian is irrelevant. - `colorvec_parabolic`: Expected to come from [SparseMatrixColorings.jl](https://github.com/gdalle/SparseMatrixColorings.jl). Allows for even faster Jacobian computation. Not necessarily required when `jac_prototype_parabolic` is given. """ From f2499a935a4e0b0315769339eeafabe357ea87fb Mon Sep 17 00:00:00 2001 From: Steve Date: Sun, 14 Sep 2025 20:06:32 -0500 Subject: [PATCH 05/22] Navierstokes hyperbolic parabolic convergence test --- ...es_convergence_implicit_sparse_jacobian.jl | 294 ++++++++++++++++++ test/test_parabolic_2d.jl | 26 ++ 2 files changed, 320 insertions(+) create mode 100644 examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl new file mode 100644 index 00000000000..9467d8e80ad --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl @@ -0,0 +1,294 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi +using SparseConnectivityTracer # For obtaining the Jacobian sparsity pattern +using SparseMatrixColorings # For obtaining the coloring vector +using OrdinaryDiffEqSDIRK, ADTypes + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations2D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), + Prandtl = prandtl_number(), + gradient_variables = GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux + +# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of +# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. +# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. +# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`. +# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`. +# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the +# `StepsizeCallback` (CFL-Condition) and less diffusion. +solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive), + volume_integral = VolumeIntegralWeakForm()) + +coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) +coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) + +# Create a uniformly refined mesh with periodic boundaries +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + periodicity = (true, false), + n_cells_max = 30_000) # set maximum capacity of tree data structure + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) +# and by the initial condition (which passes in `CompressibleEulerEquations2D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Amplitude and shift + RealT = eltype(x) + A = 0.5f0 + c = 2 + + # convenience values for trig. functions + pi_x = convert(RealT, pi) * x[1] + pi_y = convert(RealT, pi) * x[2] + pi_t = convert(RealT, pi) * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2) * (1 - exp(-A * (x[2] - 1))) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + RealT = eltype(x) + y = x[2] + + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Same settings as in `initial_condition` + # Amplitude and shift + A = 0.5f0 + c = 2 + + # convenience values for trig. functions + pi_x = convert(RealT, pi) * x[1] + pi_y = convert(RealT, pi) * x[2] + pi_t = convert(RealT, pi) * t + + # compute the manufactured solution and all necessary derivatives + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_t = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_y) * sin(pi_t) + rho_x = convert(RealT, pi) * A * cos(pi_x) * cos(pi_y) * cos(pi_t) + rho_y = -convert(RealT, pi) * A * sin(pi_x) * sin(pi_y) * cos(pi_t) + rho_xx = -convert(RealT, pi)^2 * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + rho_yy = -convert(RealT, pi)^2 * A * sin(pi_x) * cos(pi_y) * cos(pi_t) + + v1 = sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * cos(pi_t) + v1_t = -convert(RealT, pi) * sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * + sin(pi_t) + v1_x = convert(RealT, pi) * cos(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * cos(pi_t) + v1_y = sin(pi_x) * + (A * log(y + 2) * exp(-A * (y - 1)) + + (1 - exp(-A * (y - 1))) / (y + 2)) * cos(pi_t) + v1_xx = -convert(RealT, pi)^2 * sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * + cos(pi_t) + v1_xy = convert(RealT, pi) * cos(pi_x) * + (A * log(y + 2) * exp(-A * (y - 1)) + + (1 - exp(-A * (y - 1))) / (y + 2)) * cos(pi_t) + v1_yy = (sin(pi_x) * + (2 * A * exp(-A * (y - 1)) / (y + 2) - + A * A * log(y + 2) * exp(-A * (y - 1)) - + (1 - exp(-A * (y - 1))) / ((y + 2) * (y + 2))) * cos(pi_t)) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_xx = v1_xx + v2_xy = v1_xy + v2_yy = v1_yy + + p = rho * rho + p_t = 2 * rho * rho_t + p_x = 2 * rho * rho_x + p_y = 2 * rho * rho_y + p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x + p_yy = 2 * rho * rho_yy + 2 * rho_y * rho_y + + # Note this simplifies slightly because the ansatz assumes that v1 = v2 + E = p * inv_gamma_minus_one + 0.5f0 * rho * (v1^2 + v2^2) + E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2 * rho * v1 * v1_y + + # Some convenience constants + T_const = equations.gamma * inv_gamma_minus_one / Pr + inv_rho_cubed = 1 / (rho^3) + + # compute the source terms + # density equation + du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y + + # x-momentum equation + du2 = (rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y - + # stress tensor from x-direction + RealT(4) / 3 * v1_xx * mu_ + + RealT(2) / 3 * v2_xy * mu_ - + v1_yy * mu_ - + v2_xy * mu_) + # y-momentum equation + du3 = (rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2 * rho * v2 * v2_y - + # stress tensor from y-direction + v1_xy * mu_ - + v2_xx * mu_ - + RealT(4) / 3 * v2_yy * mu_ + + RealT(2) / 3 * v1_xy * mu_) + # total energy equation + du4 = (E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) - + # stress tensor and temperature gradient terms from x-direction + RealT(4) / 3 * v1_xx * v1 * mu_ + + RealT(2) / 3 * v2_xy * v1 * mu_ - + RealT(4) / 3 * v1_x * v1_x * mu_ + + RealT(2) / 3 * v2_y * v1_x * mu_ - + v1_xy * v2 * mu_ - + v2_xx * v2 * mu_ - + v1_y * v2_x * mu_ - + v2_x * v2_x * mu_ - + T_const * inv_rho_cubed * + (p_xx * rho * rho - + 2 * p_x * rho * rho_x + + 2 * p * rho_x * rho_x - + p * rho * rho_xx) * mu_ - + # stress tensor and temperature gradient terms from y-direction + v1_yy * v1 * mu_ - + v2_xy * v1 * mu_ - + v1_y * v1_y * mu_ - + v2_x * v1_y * mu_ - + RealT(4) / 3 * v2_yy * v2 * mu_ + + RealT(2) / 3 * v1_xy * v2 * mu_ - + RealT(4) / 3 * v2_y * v2_y * mu_ + + RealT(2) / 3 * v1_x * v2_y * mu_ - + T_const * inv_rho_cubed * + (p_yy * rho * rho - + 2 * p_y * rho * rho_y + + 2 * p * rho_y * rho_y - + p * rho * rho_yy) * mu_) + + return SVector(du1, du2, du3, du4) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip() do x, t, equations_parabolic + u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations_parabolic) + # This may be an `SVector` or simply a `Tuple` + return (u_cons[2] / u_cons[1], u_cons[3] / u_cons[1]) +end +heat_bc_top_bottom = Adiabatic((x, t, equations_parabolic) -> zero(eltype(x))) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, + heat_bc_top_bottom) + +# define inviscid boundary conditions +boundary_conditions = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +# define viscous boundary conditions +boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_top_bottom, + y_pos = boundary_condition_top_bottom) + +############################################################################### +### semidiscretization for sparsity detection ### + +jac_detector = TracerSparsityDetector() +# We need to construct the semidiscretization with the correct +# sparsity-detection ready datatype, which is retrieved here +jac_eltype = jacobian_eltype(real(solver), jac_detector) + +semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, + (equations, equations_parabolic), + initial_condition, solver, + uEltype = jac_eltype; + solver_parabolic = ViscousFormulationBassiRebay1(), + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic)) + +tspan = (0.0, 0.5) # Re-used for wrapping `rhs_parabolic!` below + +# Call `semidiscretize` to create the ODE problem to have access to the +# initial condition based on which the sparsity pattern is computed +ode_jac_type = semidiscretize(semi_jac_type, tspan) +u0_ode = ode_jac_type.u0 +du_ode = similar(u0_ode) + +############################################################################### +### Compute the Jacobian sparsity pattern ### + +# Wrap the `Trixi.rhs!` function to match the signature `f!(du, u)`, see +# https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity +rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type, tspan[1]) + +jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) + +# For most efficient solving we also want the coloring vector + +coloring_prob = ColoringProblem(; structure = :nonsymmetric, partition = :column) +coloring_alg = GreedyColoringAlgorithm(; decompression = :direct) +coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) +coloring_vec_parabolic = column_colors(coloring_result) + + +############################################################################### +### sparsity-aware semidiscretization and ODE ### + +# Must be called 'semi' in order for the convergence test to run successfully +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + boundary_conditions = (boundary_conditions, + boundary_conditions_parabolic), + source_terms = source_terms_navier_stokes_convergence_test) + +# Supply Jacobian prototype and coloring vector to the semidiscretization +ode_jac_sparse = semidiscretize(semi, tspan, + jac_prototype_parabolic=jac_prototype_parabolic, + colorvec_parabolic=coloring_vec_parabolic) +# using "dense" `ode = semidiscretize(semi, tspan)` +# is infeasible for larger refinement levels + +############################################################################### +# callbacks + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval = 10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode_jac_sparse, + # Default `AutoForwardDiff()` is not yet working, see + # https://github.com/trixi-framework/Trixi.jl/issues/2369 + Kvaerno4(; autodiff = AutoFiniteDiff()); + abstol = time_int_tol, reltol = time_int_tol, dt = 1e-5, + save_everystep = false, + ode_default_options()..., callback = callbacks) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 22d71b47362..fa80e42d08b 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -279,6 +279,32 @@ end end end +@trixi_testset "TreeMesh2D: elixir_navierstokes_convergence_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_navierstokes_convergence_implicit_sparse_jacobian.jl"), + tspan=(0.0, 0.5), + l2=[ + 0.0032371766320782765, + 0.006148721228526513, + 0.00407449264954935, + 0.00832496237207414 + ], + linf=[ + 0.017011038728151018, + 0.05999349864161193, + 0.01523347198959316, + 0.03580225621058908 + ]) + # 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, t)) < 1000 + end +end + @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), From 577871fea8db59267be947ed402ced90fb16ac9a Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 16 Sep 2025 10:08:02 +0200 Subject: [PATCH 06/22] Adding CI test for advection_diffusion elixir --- test/test_parabolic_2d.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index fa80e42d08b..66eddee14fc 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -188,6 +188,22 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "TreeMesh2D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", + "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), + initial_refinement_level=2, tspan=(0.0, 1.5), + l2=[0.00036653919560810473], + linf=[0.0017301065565563656]) + # 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_type, t)) < 1000 + end +end + @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"), From 2847ebe58bbe00777c108dd80516e0facba1a178 Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 16 Sep 2025 10:09:36 +0200 Subject: [PATCH 07/22] Revert "Adding CI test for advection_diffusion elixir" This reverts commit 5c4198389c5944ba473bfcde9e2446257bbd97a1. --- test/test_parabolic_2d.jl | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 66eddee14fc..fa80e42d08b 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -188,22 +188,6 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end -@trixi_testset "TreeMesh2D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", - "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), - initial_refinement_level=2, tspan=(0.0, 1.5), - l2=[0.00036653919560810473], - linf=[0.0017301065565563656]) - # 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_type, t)) < 1000 - end -end - @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"), From 0cd77a7dd737a565e16e75d3e68c47ecb3ecf39c Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 25 Sep 2025 10:32:48 +0200 Subject: [PATCH 08/22] Apply suggestions from code review --- ...tion_diffusion_implicit_sparse_jacobian.jl | 6 ++---- ...es_convergence_implicit_sparse_jacobian.jl | 21 +++++-------------- ...semidiscretization_hyperbolic_parabolic.jl | 4 ++-- 3 files changed, 9 insertions(+), 22 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index fd1c5c953f4..644256a1edd 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -18,7 +18,6 @@ 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, @@ -53,7 +52,6 @@ semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver, uEltype = jac_eltype; - solver_parabolic = ViscousFormulationBassiRebay1(), boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) @@ -68,7 +66,7 @@ du_ode = similar(u0_ode) ############################################################################### ### Compute the Jacobian sparsity pattern ### -# Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian +# Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian, see # https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitFunction # Wrap the `Trixi.rhs_parabolic!` function to match the signature `f!(du, u)`, see @@ -87,7 +85,7 @@ coloring_vec_parabolic = column_colors(coloring_result) ############################################################################### ### sparsity-aware semidiscretization and ODE ### -# Semidiscretization for actual simulation. `eEltype` is here retrieved from `solver` +# Semidiscretization for actual simulation. `uEltype` is here retrieved from `solver` semi_float_type = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl index 9467d8e80ad..b02ab2956f2 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl @@ -15,22 +15,11 @@ equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), Prandtl = prandtl_number(), gradient_variables = GradientVariablesPrimitive()) -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux - -# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of -# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`. -# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`. -# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`. -# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`. -# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the -# `StepsizeCallback` (CFL-Condition) and less diffusion. -solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive), - volume_integral = VolumeIntegralWeakForm()) - -coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y)) -coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) - -# Create a uniformly refined mesh with periodic boundaries +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 = 2, periodicity = (true, false), diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index 052957af6cc..a6d64676133 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -279,7 +279,7 @@ Optional keyword arguments: Specifies the sparsity structure of the parabolic function's Jacobian to enable e.g. efficient implicit time stepping. The [SplitODEProblem](https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitODEProblem) only expects the Jacobian to be defined on the first function it takes in, which is treated implicitly. This corresponds to the parabolic right-hand side in Trixi. - The hyperbolic right-hand side is treated explicitly, and therefore its Jacobian is irrelevant. + The hyperbolic right-hand side is expected to be treated explicitly, and therefore its Jacobian is irrelevant. - `colorvec_parabolic`: Expected to come from [SparseMatrixColorings.jl](https://github.com/gdalle/SparseMatrixColorings.jl). Allows for even faster Jacobian computation. Not necessarily required when `jac_prototype_parabolic` is given. """ @@ -343,7 +343,7 @@ Optional keyword arguments: Specifies the sparsity structure of the parabolic function's Jacobian to enable e.g. efficient implicit time stepping. The [SplitODEProblem](https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitODEProblem) only expects the Jacobian to be defined on the first function it takes in, which is treated implicitly. This corresponds to the parabolic right-hand side in Trixi. - The hyperbolic right-hand side is treated explicitly, and therefore its Jacobian is irrelevant. + The hyperbolic right-hand side is expected to be treated explicitly, and therefore its Jacobian is irrelevant. - `colorvec_parabolic`: Expected to come from [SparseMatrixColorings.jl](https://github.com/gdalle/SparseMatrixColorings.jl). Allows for even faster Jacobian computation. Not necessarily required when `jac_prototype_parabolic` is given. """ From e48fabb7e7832624731d32fbcfb074db204e8f53 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 25 Sep 2025 10:34:56 +0200 Subject: [PATCH 09/22] Apply suggestions from code review --- ...r_advection_diffusion_implicit_sparse_jacobian.jl | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 644256a1edd..24b9eebb148 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -37,9 +37,6 @@ function initial_condition_diffusive_convergence_test(x, t, end initial_condition = initial_condition_diffusive_convergence_test -boundary_conditions = boundary_condition_periodic -boundary_conditions_parabolic = boundary_condition_periodic - ############################################################################### ### semidiscretization for sparsity detection ### @@ -51,9 +48,7 @@ jac_eltype = jacobian_eltype(real(solver), jac_detector) semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver, - uEltype = jac_eltype; - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic)) + uEltype = jac_eltype) tspan = (0.0, 1.5) # Re-used for wrapping `rhs_parabolic!` below @@ -88,10 +83,7 @@ coloring_vec_parabolic = column_colors(coloring_result) # Semidiscretization for actual simulation. `uEltype` is here retrieved from `solver` semi_float_type = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), - initial_condition, solver; - solver_parabolic = ViscousFormulationBassiRebay1(), - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic)) + initial_condition, solver) # Supply Jacobian prototype and coloring vector to the semidiscretization ode_jac_sparse = semidiscretize(semi_float_type, tspan, From 385883f5a9df2535c606ba599defa0bb9d601772 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Thu, 25 Sep 2025 10:35:21 +0200 Subject: [PATCH 10/22] Update examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl --- .../elixir_navierstokes_convergence_implicit_sparse_jacobian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl index b02ab2956f2..ce7780a70ca 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl +++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl @@ -23,7 +23,7 @@ coordinates_max = (1.0, 1.0) mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2, periodicity = (true, false), - n_cells_max = 30_000) # set maximum capacity of tree data structure + n_cells_max = 30_000) # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` # since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) From e64507c3244bcb847e9f73c53951e09c049fff80 Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 16 Sep 2025 10:08:02 +0200 Subject: [PATCH 11/22] Adding CI test for advection_diffusion elixir --- test/test_parabolic_2d.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index fa80e42d08b..66eddee14fc 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -188,6 +188,22 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "TreeMesh2D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", + "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), + initial_refinement_level=2, tspan=(0.0, 1.5), + l2=[0.00036653919560810473], + linf=[0.0017301065565563656]) + # 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_type, t)) < 1000 + end +end + @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"), From b7ab829590b5156c3b5f8d80b65f12beba07a53a Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 16 Sep 2025 10:09:36 +0200 Subject: [PATCH 12/22] Revert "Adding CI test for advection_diffusion elixir" This reverts commit 5c4198389c5944ba473bfcde9e2446257bbd97a1. --- test/test_parabolic_2d.jl | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 66eddee14fc..fa80e42d08b 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -188,22 +188,6 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end -@trixi_testset "TreeMesh2D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", - "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), - initial_refinement_level=2, tspan=(0.0, 1.5), - l2=[0.00036653919560810473], - linf=[0.0017301065565563656]) - # 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_type, t)) < 1000 - end -end - @trixi_testset "TreeMesh2D: elixir_advection_diffusion.jl (LDG)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_advection_diffusion.jl"), From 734db2a2a4a08aade0619789ffc99bacaad40606 Mon Sep 17 00:00:00 2001 From: Steve Date: Fri, 26 Sep 2025 18:58:32 +0200 Subject: [PATCH 13/22] Convert adv_diffusion to 1D and N-S to P4est --- ...e_nonperiodic_implicit_sparse_jacobian.jl} | 33 ++++----- ...tion_diffusion_implicit_sparse_jacobian.jl | 12 ++-- test/test_parabolic_1d.jl | 15 +++++ test/test_parabolic_2d.jl | 67 +++++++------------ 4 files changed, 65 insertions(+), 62 deletions(-) rename examples/{tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl => p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl} (91%) rename examples/{tree_2d_dgsem => tree_1d_dgsem}/elixir_advection_diffusion_implicit_sparse_jacobian.jl (95%) diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl similarity index 91% rename from examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl rename to examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl index ce7780a70ca..54d0787f147 100644 --- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence_implicit_sparse_jacobian.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl @@ -20,10 +20,11 @@ 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 = 2, - periodicity = (true, false), - n_cells_max = 30_000) +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg=3, initial_refinement_level = 2, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + periodicity = (false, false)) # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` # since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) @@ -184,24 +185,26 @@ initial_condition = initial_condition_navier_stokes_convergence_test # BC types velocity_bc_top_bottom = NoSlip() do x, t, equations_parabolic u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations_parabolic) - # This may be an `SVector` or simply a `Tuple` - return (u_cons[2] / u_cons[1], u_cons[3] / u_cons[1]) + return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1]) end -heat_bc_top_bottom = Adiabatic((x, t, equations_parabolic) -> zero(eltype(x))) + +heat_bc_top_bottom = Adiabatic((x, t, equations_parabolic) -> 0.0) boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) +boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test) + # define inviscid boundary conditions -boundary_conditions = (; x_neg = boundary_condition_periodic, - x_pos = boundary_condition_periodic, - y_neg = boundary_condition_slip_wall, - y_pos = boundary_condition_slip_wall) +boundary_conditions = Dict(:x_neg => boundary_condition_left_right, + :x_pos => boundary_condition_left_right, + :y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall) # define viscous boundary conditions -boundary_conditions_parabolic = (; x_neg = boundary_condition_periodic, - x_pos = boundary_condition_periodic, - y_neg = boundary_condition_top_bottom, - y_pos = boundary_condition_top_bottom) +boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right, + :x_pos => boundary_condition_left_right, + :y_neg => boundary_condition_top_bottom, + :y_pos => boundary_condition_top_bottom) ############################################################################### ### semidiscretization for sparsity detection ### diff --git a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl similarity index 95% rename from examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl rename to examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 24b9eebb148..23a06d1e503 100644 --- a/examples/tree_2d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -6,22 +6,22 @@ using OrdinaryDiffEqSDIRK, ADTypes ############################################################################### # semidiscretization of the linear advection-diffusion equation -advection_velocity = (1.5, 1.0) -equations = LinearScalarAdvectionEquation2D(advection_velocity) +advection_velocity = 1.5 +equations = LinearScalarAdvectionEquation1D(advection_velocity) diffusivity() = 5.0e-2 -equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations) +equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations) solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) -coordinates_min = (-1.0, -1.0) -coordinates_max = (1.0, 1.0) +coordinates_min = -1.0 +coordinates_max = 1.0 mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 4, n_cells_max = 30_000) function initial_condition_diffusive_convergence_test(x, t, - equation::LinearScalarAdvectionEquation2D) + equation::LinearScalarAdvectionEquation1D) # Store translated coordinate for easy use of exact solution RealT = eltype(x) x_trans = x - equation.advection_velocity * t diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index d2cc65230c4..34ee1769ae3 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -95,6 +95,21 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "TreeMesh1D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", + "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), + initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, + l2=[0.05199461836224946], linf=[0.07350050454306145]) + # 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_type, t)) < 1000 + end +end + @trixi_testset "TreeMesh1D: elixir_navierstokes_convergence_periodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_navierstokes_convergence_periodic.jl"), diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index fa80e42d08b..54e4aa7210c 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -264,47 +264,6 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end -@trixi_testset "TreeMesh2D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin - @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", - "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), - initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, - l2=[4.3996570276332366e-5], linf=[7.656808122324943e-5]) - # 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_type, t)) < 1000 - end -end - -@trixi_testset "TreeMesh2D: elixir_navierstokes_convergence_implicit_sparse_jacobian.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_navierstokes_convergence_implicit_sparse_jacobian.jl"), - tspan=(0.0, 0.5), - l2=[ - 0.0032371766320782765, - 0.006148721228526513, - 0.00407449264954935, - 0.00832496237207414 - ], - linf=[ - 0.017011038728151018, - 0.05999349864161193, - 0.01523347198959316, - 0.03580225621058908 - ]) - # 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, t)) < 1000 - end -end - @trixi_testset "TreeMesh2D: elixir_navierstokes_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", "elixir_navierstokes_convergence.jl"), @@ -632,6 +591,32 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl" begin + @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + "elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl"), + tspan=(0.0, 0.5), + l2=[ + 0.0032371766320782765, + 0.006148721228526513, + 0.00407449264954935, + 0.00832496237207414 + ], + linf=[ + 0.017011038728151018, + 0.05999349864161193, + 0.01523347198959316, + 0.03580225621058908 + ]) + # 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, t)) < 1000 + end +end + @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"), From bc5624a1b57bc7356326ac3e9e44cab82bae7545 Mon Sep 17 00:00:00 2001 From: Steve Date: Mon, 29 Sep 2025 13:26:10 +0200 Subject: [PATCH 14/22] Use IMEX methods --- ...es_convergence_nonperiodic_implicit_sparse_jacobian.jl | 8 ++++---- ...elixir_advection_diffusion_implicit_sparse_jacobian.jl | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl index 54d0787f147..5a8d157de86 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl @@ -2,7 +2,7 @@ using OrdinaryDiffEqLowStorageRK using Trixi using SparseConnectivityTracer # For obtaining the Jacobian sparsity pattern using SparseMatrixColorings # For obtaining the coloring vector -using OrdinaryDiffEqSDIRK, ADTypes +using OrdinaryDiffEqBDF, ADTypes ############################################################################### # semidiscretization of the ideal compressible Navier-Stokes equations @@ -222,7 +222,7 @@ semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, boundary_conditions = (boundary_conditions, boundary_conditions_parabolic)) -tspan = (0.0, 0.5) # Re-used for wrapping `rhs_parabolic!` below +tspan = (0.0, 0.05) # Re-used for wrapping `rhs_parabolic!` below # Call `semidiscretize` to create the ODE problem to have access to the # initial condition based on which the sparsity pattern is computed @@ -280,7 +280,7 @@ time_int_tol = 1e-8 sol = solve(ode_jac_sparse, # Default `AutoForwardDiff()` is not yet working, see # https://github.com/trixi-framework/Trixi.jl/issues/2369 - Kvaerno4(; autodiff = AutoFiniteDiff()); - abstol = time_int_tol, reltol = time_int_tol, dt = 1e-5, + SBDF2(; autodiff = AutoFiniteDiff()); + abstol = time_int_tol, reltol = time_int_tol, dt = 0.0005, save_everystep = false, ode_default_options()..., callback = callbacks) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 23a06d1e503..85ba0cec1e6 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -1,7 +1,7 @@ using Trixi using SparseConnectivityTracer # For obtaining the Jacobian sparsity pattern using SparseMatrixColorings # For obtaining the coloring vector -using OrdinaryDiffEqSDIRK, ADTypes +using OrdinaryDiffEqBDF, ADTypes ############################################################################### # semidiscretization of the linear advection-diffusion equation @@ -113,7 +113,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav time_int_tol = 1.0e-9 time_abs_tol = 1.0e-9 -sol = solve(ode_jac_sparse, TRBDF2(; autodiff = AutoFiniteDiff()); - dt = 0.1, save_everystep = false, +sol = solve(ode_jac_sparse, SBDF2(; autodiff = AutoFiniteDiff()); + dt = 0.01, save_everystep = false, abstol = time_abs_tol, reltol = time_int_tol, ode_default_options()..., callback = callbacks) \ No newline at end of file From 82740724e83e177f686b2546e57341e59d683250 Mon Sep 17 00:00:00 2001 From: Steve Date: Mon, 29 Sep 2025 14:16:03 +0200 Subject: [PATCH 15/22] Add demonstration of adv_velocity=0 --- ...tion_diffusion_implicit_sparse_jacobian.jl | 41 ++++++++++++++++--- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 85ba0cec1e6..068d1651e6a 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -7,9 +7,9 @@ using OrdinaryDiffEqBDF, ADTypes # semidiscretization of the linear advection-diffusion equation advection_velocity = 1.5 -equations = LinearScalarAdvectionEquation1D(advection_velocity) +equations_hyperbolic = LinearScalarAdvectionEquation1D(advection_velocity) diffusivity() = 5.0e-2 -equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations) +equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations_hyperbolic) solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) @@ -46,7 +46,7 @@ jac_detector = TracerSparsityDetector() jac_eltype = jacobian_eltype(real(solver), jac_detector) semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, - (equations, equations_parabolic), + (equations_hyperbolic, equations_parabolic), initial_condition, solver, uEltype = jac_eltype) @@ -63,6 +63,9 @@ du_ode = similar(u0_ode) # Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian, see # https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitFunction +# Although we do sparsity detection on the entire RHS (since semi_jac_type depends on both equations and +# equations_parabolic), this is equivalent to doing sparsity detection on the diffusion problem alone +# (see next section for a demonstration of this) # Wrap the `Trixi.rhs_parabolic!` function to match the signature `f!(du, u)`, see # https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity @@ -77,12 +80,40 @@ coloring_alg = GreedyColoringAlgorithm(; decompression = :direct) coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) coloring_vec_parabolic = column_colors(coloring_result) +############################################################################### +### Compare sparsity pattern based on hyperbolic-parabolic ### +### problem with sparsity pattern of parabolic-only problem ### + +# We construct a semidiscretization just as we did previously, but this time with advection_velocity=0 +equations_hyperbolic_no_advection = LinearScalarAdvectionEquation1D(0) +equations_parabolic_no_advection = LaplaceDiffusion1D(diffusivity(), equations_hyperbolic_no_advection) + +semi_jac_type_no_advection = SemidiscretizationHyperbolicParabolic(mesh, + (equations_hyperbolic_no_advection, + equations_parabolic_no_advection), + initial_condition, solver, + uEltype = jac_eltype) +ode_jac_type_no_advection = semidiscretize(semi_jac_type_no_advection, tspan) + +# Do sparsity detection on our semidiscretization with advection turned off +rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type_no_advection, tspan[1]) + +jac_prototype_parabolic_no_advection = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) + +coloring_result_no_advection = coloring(jac_prototype_parabolic_no_advection, coloring_prob, coloring_alg) +coloring_vec_parabolic_no_advection = column_colors(coloring_result) + +# Given that the stencil for parabolic solvers are always larger than those of hyperbolic solvers, +# the sparsity detection will never depend on the hyperbolic part of the problem +@assert(jac_prototype_parabolic == jac_prototype_parabolic_no_advection) +@assert(coloring_vec_parabolic == coloring_vec_parabolic_no_advection) + ############################################################################### ### sparsity-aware semidiscretization and ODE ### # Semidiscretization for actual simulation. `uEltype` is here retrieved from `solver` semi_float_type = SemidiscretizationHyperbolicParabolic(mesh, - (equations, equations_parabolic), + (equations_hyperbolic, equations_parabolic), initial_condition, solver) # Supply Jacobian prototype and coloring vector to the semidiscretization @@ -116,4 +147,4 @@ time_abs_tol = 1.0e-9 sol = solve(ode_jac_sparse, SBDF2(; autodiff = AutoFiniteDiff()); dt = 0.01, save_everystep = false, abstol = time_abs_tol, reltol = time_int_tol, - ode_default_options()..., callback = callbacks) \ No newline at end of file + ode_default_options()..., callback = callbacks) From 89ce298d8242c2814b45274053a8736ad6db4d41 Mon Sep 17 00:00:00 2001 From: Steve Date: Mon, 29 Sep 2025 14:26:34 +0200 Subject: [PATCH 16/22] Fix tests --- test/test_parabolic_1d.jl | 2 +- test/test_parabolic_2d.jl | 27 +++++++++++---------------- 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index 34ee1769ae3..fc05b460022 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -99,7 +99,7 @@ end @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, - l2=[0.05199461836224946], linf=[0.07350050454306145]) + l2=[0.05240130204342638], linf=[0.07407444680136666]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 54e4aa7210c..944c13020f9 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -592,29 +592,24 @@ end end @trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl" begin - @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", + @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl"), - tspan=(0.0, 0.5), + tspan=(0.0, 0.003), l2=[ - 0.0032371766320782765, - 0.006148721228526513, - 0.00407449264954935, - 0.00832496237207414 + 3.1541556275071338e-6, + 8.184640302981804e-6, + 8.062391655601536e-6, + 4.7582512486365587e-5 ], linf=[ - 0.017011038728151018, - 0.05999349864161193, - 0.01523347198959316, - 0.03580225621058908 + 1.5362999041812486e-5, + 3.765373955166851e-5, + 4.380899179601272e-5, + 0.000278756949764869 ]) # 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, t)) < 1000 - end + @test_allocations(Trixi.rhs!, semi, sol, 1000) end @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin From 57d7d978092c2d613b99d249894ae87195f61a87 Mon Sep 17 00:00:00 2001 From: Daniel Doehring Date: Mon, 29 Sep 2025 16:21:09 +0200 Subject: [PATCH 17/22] Apply suggestions from code review --- .../semidiscretization_hyperbolic_parabolic.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl index a6d64676133..7943cbc60ea 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl @@ -316,10 +316,6 @@ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan; else # We could also construct an `ODEFunction` without the Jacobian here, # but we stick to the more light-weight direct in-place function `rhs_parabolic!`. - - # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the - # first function implicitly and the second one explicitly. Thus, we pass the - # stiffer parabolic function first. return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) end end @@ -381,10 +377,6 @@ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan, else # We could also construct an `ODEFunction` without the Jacobian here, # but we stick to the more light-weight direct in-place function `rhs_parabolic!`. - - # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the - # first function implicitly and the second one explicitly. Thus, we pass the - # stiffer parabolic function first. return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi) end end From a2586114306266fa8948d4476fb568c5f267d8b0 Mon Sep 17 00:00:00 2001 From: Steve Date: Mon, 29 Sep 2025 20:45:26 +0200 Subject: [PATCH 18/22] PR comments --- ...tion_diffusion_implicit_sparse_jacobian.jl | 50 ------------ test/test_parabolic_2d.jl | 17 ++-- test/test_unit.jl | 81 +++++++++++++++++++ 3 files changed, 90 insertions(+), 58 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 068d1651e6a..1c0a52fda4c 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -58,56 +58,6 @@ ode_jac_type = semidiscretize(semi_jac_type, tspan) u0_ode = ode_jac_type.u0 du_ode = similar(u0_ode) -############################################################################### -### Compute the Jacobian sparsity pattern ### - -# Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian, see -# https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitFunction -# Although we do sparsity detection on the entire RHS (since semi_jac_type depends on both equations and -# equations_parabolic), this is equivalent to doing sparsity detection on the diffusion problem alone -# (see next section for a demonstration of this) - -# Wrap the `Trixi.rhs_parabolic!` function to match the signature `f!(du, u)`, see -# https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity -rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type, tspan[1]) - -jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) - -# For most efficient solving we also want the coloring vector - -coloring_prob = ColoringProblem(; structure = :nonsymmetric, partition = :column) -coloring_alg = GreedyColoringAlgorithm(; decompression = :direct) -coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) -coloring_vec_parabolic = column_colors(coloring_result) - -############################################################################### -### Compare sparsity pattern based on hyperbolic-parabolic ### -### problem with sparsity pattern of parabolic-only problem ### - -# We construct a semidiscretization just as we did previously, but this time with advection_velocity=0 -equations_hyperbolic_no_advection = LinearScalarAdvectionEquation1D(0) -equations_parabolic_no_advection = LaplaceDiffusion1D(diffusivity(), equations_hyperbolic_no_advection) - -semi_jac_type_no_advection = SemidiscretizationHyperbolicParabolic(mesh, - (equations_hyperbolic_no_advection, - equations_parabolic_no_advection), - initial_condition, solver, - uEltype = jac_eltype) -ode_jac_type_no_advection = semidiscretize(semi_jac_type_no_advection, tspan) - -# Do sparsity detection on our semidiscretization with advection turned off -rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type_no_advection, tspan[1]) - -jac_prototype_parabolic_no_advection = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) - -coloring_result_no_advection = coloring(jac_prototype_parabolic_no_advection, coloring_prob, coloring_alg) -coloring_vec_parabolic_no_advection = column_colors(coloring_result) - -# Given that the stencil for parabolic solvers are always larger than those of hyperbolic solvers, -# the sparsity detection will never depend on the hyperbolic part of the problem -@assert(jac_prototype_parabolic == jac_prototype_parabolic_no_advection) -@assert(coloring_vec_parabolic == coloring_vec_parabolic_no_advection) - ############################################################################### ### sparsity-aware semidiscretization and ODE ### diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 944c13020f9..edb0d6d3f26 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -595,17 +595,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl"), tspan=(0.0, 0.003), + initial_refinement_level=1, l2=[ - 3.1541556275071338e-6, - 8.184640302981804e-6, - 8.062391655601536e-6, - 4.7582512486365587e-5 + 4.588219785453587e-5, + 0.00011485891993215742, + 0.00011457767113486697, + 0.0007063058724082755 ], linf=[ - 1.5362999041812486e-5, - 3.765373955166851e-5, - 4.380899179601272e-5, - 0.000278756949764869 + 0.00017051082943986273, + 0.0005007782853971854, + 0.0005309051938069409, + 0.0036206287734668052 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unit.jl b/test/test_unit.jl index 54403a3e3c2..33a4a76ad46 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -2692,6 +2692,87 @@ end @test isapprox(jac_finite_diff, Matrix(jac_sparse_finite_diff); rtol = 5e-8) @test isapprox(sparse(jac_finite_diff), jac_sparse_finite_diff; rtol = 5e-8) end + +@testset "Parabolic-Hyperbolic Problem Sparsity Pattern " begin + ############################################################################### + ### equations, solver, mesh ### + + advection_velocity = 1.5 + equations_hyperbolic = LinearScalarAdvectionEquation1D(advection_velocity) + diffusivity() = 5.0e-2 + equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations_hyperbolic) + + solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) + + coordinates_min = -1.0 + coordinates_max = 1.0 + + mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 30_000) + + ############################################################################### + ### semidiscretization for sparsity detection ### + + jac_detector = TracerSparsityDetector() + # We need to construct the semidiscretization with the correct + # sparsity-detection ready datatype, which is retrieved here + jac_eltype = jacobian_eltype(real(solver), jac_detector) + + # Semidiscretization for sparsity pattern detection + semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, (equations_hyperbolic, equations_parabolic), + initial_condition_convergence_test, + solver, + uEltype = jac_eltype) # Need to supply Jacobian element type + + tspan = (0.0, 1.5) # Re-used for wrapping `rhs` below + + # Call `semidiscretize` to create the ODE problem to have access to the + # initial condition based on which the sparsity pattern is computed + ode_jac_type = semidiscretize(semi_jac_type, tspan) + u0_ode = ode_jac_type.u0 + du_ode = similar(u0_ode) + + ############################################################################### + ### Compute the Jacobian sparsity pattern ### + + # Only the parabolic part of the `SplitODEProblem` is treated implicitly so we only need the parabolic Jacobian, see + # https://docs.sciml.ai/DiffEqDocs/stable/types/split_ode_types/#SciMLBase.SplitFunction + # Although we do sparsity detection on the entire RHS (since semi_jac_type depends on both equations and + # equations_parabolic), this is equivalent to doing sparsity detection on the diffusion problem alone + # (see next section for a demonstration of this) + + # Wrap the `Trixi.rhs_parabolic!` function to match the signature `f!(du, u)`, see + # https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity + rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type, tspan[1]) + + jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) + ############################################################################### + ### Compare sparsity pattern based on hyperbolic-parabolic ### + ### problem with sparsity pattern of parabolic-only problem ### + + # We construct a semidiscretization just as we did previously, but this time with advection_velocity=0 + equations_hyperbolic_no_advection = LinearScalarAdvectionEquation1D(0) + equations_parabolic_no_advection = LaplaceDiffusion1D(diffusivity(), equations_hyperbolic_no_advection) + + semi_jac_type_no_advection = SemidiscretizationHyperbolicParabolic(mesh, + (equations_hyperbolic_no_advection, + equations_parabolic_no_advection), + initial_condition, solver, + uEltype = jac_eltype) + ode_jac_type_no_advection = semidiscretize(semi_jac_type_no_advection, tspan) + + # Do sparsity detection on our semidiscretization with advection turned off + rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type_no_advection, tspan[1]) + + jac_prototype_parabolic_no_advection = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) + + # Given that the stencil for parabolic solvers are always larger than those of hyperbolic solvers, + # the sparsity detection will never depend on the hyperbolic part of the problem + @assert(jac_prototype_parabolic == jac_prototype_parabolic_no_advection) + + @test isapprox(jac_prototype_parabolic, jac_prototype_parabolic_no_advection; rtol = 5e-8) +end end end #module From 2fdc31578e84dc5c2667face191d31d1f0b8d6a8 Mon Sep 17 00:00:00 2001 From: Steve Date: Thu, 2 Oct 2025 16:28:12 +0200 Subject: [PATCH 19/22] PR comments2 --- Project.toml | 9 ++++++++- ...gence_nonperiodic_implicit_sparse_jacobian.jl | 2 +- ...vection_diffusion_implicit_sparse_jacobian.jl | 16 ++++++++++++++++ test/test_parabolic_1d.jl | 2 +- test/test_unit.jl | 6 ++---- 5 files changed, 28 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index a7a0b3fc48f..807f2bdbfe4 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Michael Schlottke-Lakemper ", " version = "0.13.9-DEV" [deps] +ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" @@ -15,16 +16,21 @@ DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" +OrdinaryDiffEqIMEXMultistep = "9f002381-b378-40b7-97a6-27a27c83f129" +OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" @@ -37,6 +43,8 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" +SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" @@ -59,7 +67,6 @@ Convex = "f65535da-76fb-5f13-bab9-19810c17039a" ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" [extensions] TrixiCUDAExt = "CUDA" diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl index 5a8d157de86..a148a809d61 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl @@ -24,7 +24,7 @@ trees_per_dimension = (4, 4) mesh = P4estMesh(trees_per_dimension, polydeg=3, initial_refinement_level = 2, coordinates_min=coordinates_min, coordinates_max=coordinates_max, - periodicity = (false, false)) + periodicity = false) # Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` # since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) diff --git a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl index 1c0a52fda4c..e098f1484c2 100644 --- a/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl +++ b/examples/tree_1d_dgsem/elixir_advection_diffusion_implicit_sparse_jacobian.jl @@ -58,6 +58,22 @@ ode_jac_type = semidiscretize(semi_jac_type, tspan) u0_ode = ode_jac_type.u0 du_ode = similar(u0_ode) +############################################################################### +### Compute the Jacobian sparsity pattern ### + +# Wrap the `Trixi.rhs!` function to match the signature `f!(du, u)`, see +# https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity +rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type, tspan[1]) + +jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) + +# For most efficient solving we also want the coloring vector + +coloring_prob = ColoringProblem(; structure = :nonsymmetric, partition = :column) +coloring_alg = GreedyColoringAlgorithm(; decompression = :direct) +coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) +coloring_vec_parabolic = column_colors(coloring_result) + ############################################################################### ### sparsity-aware semidiscretization and ODE ### diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl index fc05b460022..6058880f9bf 100644 --- a/test/test_parabolic_1d.jl +++ b/test/test_parabolic_1d.jl @@ -98,7 +98,7 @@ end @trixi_testset "TreeMesh1D: elixir_advection_diffusion_implicit_sparse_jacobian.jl" begin @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion_implicit_sparse_jacobian.jl"), - initial_refinement_level=4, tspan=(0.0, 0.4), polydeg=3, + tspan=(0.0, 0.4), l2=[0.05240130204342638], linf=[0.07407444680136666]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unit.jl b/test/test_unit.jl index 33a4a76ad46..3c5f2ec33fa 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -2758,7 +2758,7 @@ end semi_jac_type_no_advection = SemidiscretizationHyperbolicParabolic(mesh, (equations_hyperbolic_no_advection, equations_parabolic_no_advection), - initial_condition, solver, + initial_condition_convergence_test, solver, uEltype = jac_eltype) ode_jac_type_no_advection = semidiscretize(semi_jac_type_no_advection, tspan) @@ -2769,9 +2769,7 @@ end # Given that the stencil for parabolic solvers are always larger than those of hyperbolic solvers, # the sparsity detection will never depend on the hyperbolic part of the problem - @assert(jac_prototype_parabolic == jac_prototype_parabolic_no_advection) - - @test isapprox(jac_prototype_parabolic, jac_prototype_parabolic_no_advection; rtol = 5e-8) + @test jac_prototype_parabolic == jac_prototype_parabolic_no_advection end end From 6ddfdad56ac685381d847aff0040da4ef95aa88b Mon Sep 17 00:00:00 2001 From: Steve Date: Thu, 2 Oct 2025 17:00:15 +0200 Subject: [PATCH 20/22] Undo overwrite of Project.toml --- Project.toml | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/Project.toml b/Project.toml index 807f2bdbfe4..a7a0b3fc48f 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ authors = ["Michael Schlottke-Lakemper ", " version = "0.13.9-DEV" [deps] -ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" @@ -16,21 +15,16 @@ DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" -FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e" -LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Octavian = "6fd5a793-0b7e-452c-907f-f8bfe9c57db4" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" -OrdinaryDiffEqIMEXMultistep = "9f002381-b378-40b7-97a6-27a27c83f129" -OrdinaryDiffEqLowStorageRK = "b0944070-b475-4768-8dec-fb6eb410534d" P4est = "7d669430-f675-4ae7-b43e-fab78ec5a902" Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" @@ -43,8 +37,6 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" -SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" @@ -67,6 +59,7 @@ Convex = "f65535da-76fb-5f13-bab9-19810c17039a" ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" [extensions] TrixiCUDAExt = "CUDA" From 23e6c17e5b4e042b3ae33652d6d9a5d5891fcf29 Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 14 Oct 2025 18:31:08 +0200 Subject: [PATCH 21/22] Remove the navierstokes elixir and test --- ...ce_nonperiodic_implicit_sparse_jacobian.jl | 286 ------------------ test/test_parabolic_2d.jl | 22 -- 2 files changed, 308 deletions(-) delete mode 100644 examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl deleted file mode 100644 index a148a809d61..00000000000 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl +++ /dev/null @@ -1,286 +0,0 @@ -using OrdinaryDiffEqLowStorageRK -using Trixi -using SparseConnectivityTracer # For obtaining the Jacobian sparsity pattern -using SparseMatrixColorings # For obtaining the coloring vector -using OrdinaryDiffEqBDF, ADTypes - -############################################################################### -# semidiscretization of the ideal compressible Navier-Stokes equations - -prandtl_number() = 0.72 -mu() = 0.01 - -equations = CompressibleEulerEquations2D(1.4) -equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(), - Prandtl = prandtl_number(), - gradient_variables = GradientVariablesPrimitive()) - -solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs) - -coordinates_min = (-1.0, -1.0) -coordinates_max = (1.0, 1.0) - -trees_per_dimension = (4, 4) -mesh = P4estMesh(trees_per_dimension, - polydeg=3, initial_refinement_level = 2, - coordinates_min=coordinates_min, coordinates_max=coordinates_max, - periodicity = false) - -# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D` -# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`) -# and by the initial condition (which passes in `CompressibleEulerEquations2D`). -# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) -function initial_condition_navier_stokes_convergence_test(x, t, equations) - # Amplitude and shift - RealT = eltype(x) - A = 0.5f0 - c = 2 - - # convenience values for trig. functions - pi_x = convert(RealT, pi) * x[1] - pi_y = convert(RealT, pi) * x[2] - pi_t = convert(RealT, pi) * t - - rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - v1 = sin(pi_x) * log(x[2] + 2) * (1 - exp(-A * (x[2] - 1))) * cos(pi_t) - v2 = v1 - p = rho^2 - - return prim2cons(SVector(rho, v1, v2, p), equations) -end - -@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) - RealT = eltype(x) - y = x[2] - - # TODO: parabolic - # we currently need to hardcode these parameters until we fix the "combined equation" issue - # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 - inv_gamma_minus_one = inv(equations.gamma - 1) - Pr = prandtl_number() - mu_ = mu() - - # Same settings as in `initial_condition` - # Amplitude and shift - A = 0.5f0 - c = 2 - - # convenience values for trig. functions - pi_x = convert(RealT, pi) * x[1] - pi_y = convert(RealT, pi) * x[2] - pi_t = convert(RealT, pi) * t - - # compute the manufactured solution and all necessary derivatives - rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_t = -convert(RealT, pi) * A * sin(pi_x) * cos(pi_y) * sin(pi_t) - rho_x = convert(RealT, pi) * A * cos(pi_x) * cos(pi_y) * cos(pi_t) - rho_y = -convert(RealT, pi) * A * sin(pi_x) * sin(pi_y) * cos(pi_t) - rho_xx = -convert(RealT, pi)^2 * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - rho_yy = -convert(RealT, pi)^2 * A * sin(pi_x) * cos(pi_y) * cos(pi_t) - - v1 = sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * cos(pi_t) - v1_t = -convert(RealT, pi) * sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * - sin(pi_t) - v1_x = convert(RealT, pi) * cos(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * cos(pi_t) - v1_y = sin(pi_x) * - (A * log(y + 2) * exp(-A * (y - 1)) + - (1 - exp(-A * (y - 1))) / (y + 2)) * cos(pi_t) - v1_xx = -convert(RealT, pi)^2 * sin(pi_x) * log(y + 2) * (1 - exp(-A * (y - 1))) * - cos(pi_t) - v1_xy = convert(RealT, pi) * cos(pi_x) * - (A * log(y + 2) * exp(-A * (y - 1)) + - (1 - exp(-A * (y - 1))) / (y + 2)) * cos(pi_t) - v1_yy = (sin(pi_x) * - (2 * A * exp(-A * (y - 1)) / (y + 2) - - A * A * log(y + 2) * exp(-A * (y - 1)) - - (1 - exp(-A * (y - 1))) / ((y + 2) * (y + 2))) * cos(pi_t)) - v2 = v1 - v2_t = v1_t - v2_x = v1_x - v2_y = v1_y - v2_xx = v1_xx - v2_xy = v1_xy - v2_yy = v1_yy - - p = rho * rho - p_t = 2 * rho * rho_t - p_x = 2 * rho * rho_x - p_y = 2 * rho * rho_y - p_xx = 2 * rho * rho_xx + 2 * rho_x * rho_x - p_yy = 2 * rho * rho_yy + 2 * rho_y * rho_y - - # Note this simplifies slightly because the ansatz assumes that v1 = v2 - E = p * inv_gamma_minus_one + 0.5f0 * rho * (v1^2 + v2^2) - E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2 * rho * v1 * v1_t - E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2 * rho * v1 * v1_x - E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2 * rho * v1 * v1_y - - # Some convenience constants - T_const = equations.gamma * inv_gamma_minus_one / Pr - inv_rho_cubed = 1 / (rho^3) - - # compute the source terms - # density equation - du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y - - # x-momentum equation - du2 = (rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 - + 2 * rho * v1 * v1_x - + rho_y * v1 * v2 - + rho * v1_y * v2 - + rho * v1 * v2_y - - # stress tensor from x-direction - RealT(4) / 3 * v1_xx * mu_ + - RealT(2) / 3 * v2_xy * mu_ - - v1_yy * mu_ - - v2_xy * mu_) - # y-momentum equation - du3 = (rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 - + rho * v1_x * v2 - + rho * v1 * v2_x - + rho_y * v2^2 - + 2 * rho * v2 * v2_y - - # stress tensor from y-direction - v1_xy * mu_ - - v2_xx * mu_ - - RealT(4) / 3 * v2_yy * mu_ + - RealT(2) / 3 * v1_xy * mu_) - # total energy equation - du4 = (E_t + v1_x * (E + p) + v1 * (E_x + p_x) - + v2_y * (E + p) + v2 * (E_y + p_y) - - # stress tensor and temperature gradient terms from x-direction - RealT(4) / 3 * v1_xx * v1 * mu_ + - RealT(2) / 3 * v2_xy * v1 * mu_ - - RealT(4) / 3 * v1_x * v1_x * mu_ + - RealT(2) / 3 * v2_y * v1_x * mu_ - - v1_xy * v2 * mu_ - - v2_xx * v2 * mu_ - - v1_y * v2_x * mu_ - - v2_x * v2_x * mu_ - - T_const * inv_rho_cubed * - (p_xx * rho * rho - - 2 * p_x * rho * rho_x + - 2 * p * rho_x * rho_x - - p * rho * rho_xx) * mu_ - - # stress tensor and temperature gradient terms from y-direction - v1_yy * v1 * mu_ - - v2_xy * v1 * mu_ - - v1_y * v1_y * mu_ - - v2_x * v1_y * mu_ - - RealT(4) / 3 * v2_yy * v2 * mu_ + - RealT(2) / 3 * v1_xy * v2 * mu_ - - RealT(4) / 3 * v2_y * v2_y * mu_ + - RealT(2) / 3 * v1_x * v2_y * mu_ - - T_const * inv_rho_cubed * - (p_yy * rho * rho - - 2 * p_y * rho * rho_y + - 2 * p * rho_y * rho_y - - p * rho * rho_yy) * mu_) - - return SVector(du1, du2, du3, du4) -end - -initial_condition = initial_condition_navier_stokes_convergence_test - -# BC types -velocity_bc_top_bottom = NoSlip() do x, t, equations_parabolic - u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations_parabolic) - return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1]) -end - -heat_bc_top_bottom = Adiabatic((x, t, equations_parabolic) -> 0.0) -boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, - heat_bc_top_bottom) - -boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test) - -# define inviscid boundary conditions -boundary_conditions = Dict(:x_neg => boundary_condition_left_right, - :x_pos => boundary_condition_left_right, - :y_neg => boundary_condition_slip_wall, - :y_pos => boundary_condition_slip_wall) - -# define viscous boundary conditions -boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right, - :x_pos => boundary_condition_left_right, - :y_neg => boundary_condition_top_bottom, - :y_pos => boundary_condition_top_bottom) - -############################################################################### -### semidiscretization for sparsity detection ### - -jac_detector = TracerSparsityDetector() -# We need to construct the semidiscretization with the correct -# sparsity-detection ready datatype, which is retrieved here -jac_eltype = jacobian_eltype(real(solver), jac_detector) - -semi_jac_type = SemidiscretizationHyperbolicParabolic(mesh, - (equations, equations_parabolic), - initial_condition, solver, - uEltype = jac_eltype; - solver_parabolic = ViscousFormulationBassiRebay1(), - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic)) - -tspan = (0.0, 0.05) # Re-used for wrapping `rhs_parabolic!` below - -# Call `semidiscretize` to create the ODE problem to have access to the -# initial condition based on which the sparsity pattern is computed -ode_jac_type = semidiscretize(semi_jac_type, tspan) -u0_ode = ode_jac_type.u0 -du_ode = similar(u0_ode) - -############################################################################### -### Compute the Jacobian sparsity pattern ### - -# Wrap the `Trixi.rhs!` function to match the signature `f!(du, u)`, see -# https://adrianhill.de/SparseConnectivityTracer.jl/stable/user/api/#ADTypes.jacobian_sparsity -rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type, tspan[1]) - -jac_prototype_parabolic = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector) - -# For most efficient solving we also want the coloring vector - -coloring_prob = ColoringProblem(; structure = :nonsymmetric, partition = :column) -coloring_alg = GreedyColoringAlgorithm(; decompression = :direct) -coloring_result = coloring(jac_prototype_parabolic, coloring_prob, coloring_alg) -coloring_vec_parabolic = column_colors(coloring_result) - - -############################################################################### -### sparsity-aware semidiscretization and ODE ### - -# Must be called 'semi' in order for the convergence test to run successfully -semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), - initial_condition, solver; - boundary_conditions = (boundary_conditions, - boundary_conditions_parabolic), - source_terms = source_terms_navier_stokes_convergence_test) - -# Supply Jacobian prototype and coloring vector to the semidiscretization -ode_jac_sparse = semidiscretize(semi, tspan, - jac_prototype_parabolic=jac_prototype_parabolic, - colorvec_parabolic=coloring_vec_parabolic) -# using "dense" `ode = semidiscretize(semi, tspan)` -# is infeasible for larger refinement levels - -############################################################################### -# callbacks - -summary_callback = SummaryCallback() -alive_callback = AliveCallback(alive_interval = 10) -analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval) -callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) - -############################################################################### -# run the simulation - -time_int_tol = 1e-8 -sol = solve(ode_jac_sparse, - # Default `AutoForwardDiff()` is not yet working, see - # https://github.com/trixi-framework/Trixi.jl/issues/2369 - SBDF2(; autodiff = AutoFiniteDiff()); - abstol = time_int_tol, reltol = time_int_tol, dt = 0.0005, - save_everystep = false, - ode_default_options()..., callback = callbacks) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index edb0d6d3f26..0d23b43ef4b 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -591,28 +591,6 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end -@trixi_testset "P4estMesh2D: elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", - "elixir_navierstokes_convergence_nonperiodic_implicit_sparse_jacobian.jl"), - tspan=(0.0, 0.003), - initial_refinement_level=1, - l2=[ - 4.588219785453587e-5, - 0.00011485891993215742, - 0.00011457767113486697, - 0.0007063058724082755 - ], - linf=[ - 0.00017051082943986273, - 0.0005007782853971854, - 0.0005309051938069409, - 0.0036206287734668052 - ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - @test_allocations(Trixi.rhs!, semi, sol, 1000) -end - @trixi_testset "P4estMesh2D: elixir_navierstokes_lid_driven_cavity.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_navierstokes_lid_driven_cavity.jl"), From 4fc9196f1cf344eaccd451ebce746b8fb65c6596 Mon Sep 17 00:00:00 2001 From: Steve Date: Tue, 14 Oct 2025 18:43:04 +0200 Subject: [PATCH 22/22] Update unit test to actually include rhs! --- test/test_unit.jl | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/test/test_unit.jl b/test/test_unit.jl index 3c5f2ec33fa..9106abc2996 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -2693,7 +2693,25 @@ end @test isapprox(sparse(jac_finite_diff), jac_sparse_finite_diff; rtol = 5e-8) end -@testset "Parabolic-Hyperbolic Problem Sparsity Pattern " begin +@testset "Parabolic-Hyperbolic Problem Sparsity Pattern" begin + + function rhs_hyperbolic_parabolic!(du_ode, u_ode, + semi::SemidiscretizationHyperbolicParabolic, t) + Trixi.@trixi_timeit timer() "rhs_hyperbolic_parabolic!" begin + # Implementation of split ODE problem in OrdinaryDiffEq + du_para = similar(du_ode) # This obviously allocates + rhs!(du_ode, u_ode, semi, t) + rhs_parabolic!(du_para, u_ode, semi, t) + + Trixi.@threaded for i in eachindex(du_ode) + # Try to enable optimizations due to `muladd` by avoiding `+=` + # https://github.com/trixi-framework/Trixi.jl/pull/2480#discussion_r2224531702 + du_ode[i] = du_ode[i] + du_para[i] + end + end + return nothing + end + ############################################################################### ### equations, solver, mesh ### @@ -2763,7 +2781,7 @@ end ode_jac_type_no_advection = semidiscretize(semi_jac_type_no_advection, tspan) # Do sparsity detection on our semidiscretization with advection turned off - rhs_parabolic_wrapped! = (du_ode, u0_ode) -> Trixi.rhs_parabolic!(du_ode, u0_ode, semi_jac_type_no_advection, tspan[1]) + rhs_hyper_para! = (du_ode, u0_ode) -> rhs_hyperbolic_parabolic!(du_ode, u0_ode, semi_jac_type_no_advection, tspan[1]) jac_prototype_parabolic_no_advection = jacobian_sparsity(rhs_parabolic_wrapped!, du_ode, u0_ode, jac_detector)