Update to OrdinaryDiffEq.jl v7 and related SciML packages

Project.toml

@@ -61,6 +61,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"
+OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 
@@ -71,6 +72,7 @@ TrixiConvexECOSExt = ["Convex", "ECOS"]
 TrixiMakieExt = "Makie"
 TrixiNLsolveExt = "NLsolve"
 TrixiPlotsExt = "Plots"
+TrixiOrdinaryDiffEqCoreExt = "OrdinaryDiffEqCore"
 TrixiSparseConnectivityTracerExt = "SparseConnectivityTracer"
 
 [compat]
@@ -84,7 +86,7 @@ ConstructionBase = "1.5.8"
 Convex = "0.16"
 DataStructures = "0.18.15, 0.19"
 DelimitedFiles = "1"
-DiffEqBase = "6.184"
+DiffEqBase = "6.184, 7"
 DiffEqCallbacks = "2.35, 3, 4"
 Downloads = "1.6"
 ECOS = "1.1.2"
@@ -103,16 +105,17 @@ MuladdMacro = "0.2.4"
 NLsolve = "4.5.1"
 Octavian = "0.3.28"
 OffsetArrays = "1.13"
+OrdinaryDiffEqCore = "1.26, 2, 3, 4.2"
 P4est = "0.4.12"
 Plots = "1.38.13"
 Polyester = "=0.7.16, 0.7.18"
 PrecompileTools = "1.2.1"
 Preferences = "1.5"
 Printf = "1"
 RecipesBase = "1.3.4"
-RecursiveArrayTools = "3.37"
+RecursiveArrayTools = "3.37, 4"
 Reexport = "1.2.2"
-SciMLBase = "2.128"
+SciMLBase = "2.128, 3"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 SparseConnectivityTracer = "1.0.1"

benchmark/Project.toml

@@ -6,6 +6,6 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 
 [compat]
 BenchmarkTools = "0.5, 0.7, 1.0"
-OrdinaryDiffEq = "5.65, 6"
+OrdinaryDiffEq = "5.65, 6, 7"
 PkgBenchmark = "0.2.10"
 Trixi = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16"

docs/Project.toml

@@ -41,11 +41,11 @@ LaTeXStrings = "1.2"
 Literate = "2.9"
 Measurements = "2.5"
 NLsolve = "4.5.1"
-OrdinaryDiffEqLowOrderRK = "1.2"
-OrdinaryDiffEqLowStorageRK = "1.2"
-OrdinaryDiffEqSDIRK = "1.1"
-OrdinaryDiffEqSSPRK = "1.2"
-OrdinaryDiffEqTsit5 = "1.1"
+OrdinaryDiffEqLowOrderRK = "1.2, 2"
+OrdinaryDiffEqLowStorageRK = "1.2, 2, 3"
+OrdinaryDiffEqSDIRK = "1.1, 2"
+OrdinaryDiffEqSSPRK = "1.2, 2"
+OrdinaryDiffEqTsit5 = "1.1, 2"
 Plots = "1.9"
 SparseConnectivityTracer = "1.0.1"
 SparseMatrixColorings = "0.4.21"

examples/tree_2d_dgsem/elixir_euler_colliding_flow_amr.jl

@@ -119,5 +119,5 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6, 5.0e-
 ###############################################################################
 # run the simulation
 # use adaptive time stepping based on error estimates, time step roughly dt = 5e-3
-sol = solve(ode, SSPRK43(stage_limiter!);
+sol = solve(ode, SSPRK43(; stage_limiter!);
             ode_default_options()..., callback = callbacks);

ext/TrixiOrdinaryDiffEqCoreExt.jl

@@ -0,0 +1,113 @@
+module TrixiOrdinaryDiffEqCoreExt
+
+import Trixi: load_controller!, store_controller!
+import OrdinaryDiffEqCore: OrdinaryDiffEqCore, PIController, PIDController
+import HDF5: attributes
+
+@static if Base.pkgversion(OrdinaryDiffEqCore) >= v"4"
+    import OrdinaryDiffEqCore: PIControllerCache, PIDControllerCache
+end
+
+# Support to load controller.
+# OrdinaryDiffEqCore < v4: The state of PIController and PIDController lives on the
+# integrator (integrator.qold) rather than on the controller struct.
+# OrdinaryDiffEqCore >= v4: the legacy controllers are replaced by PIControllerCache and
+# PIDControllerCache, and integrator.qold no longer exists.
+
+# OrdinaryDiffEqCore < v4, PI controller:
+# Previous error estimate stored as integrator.qold.
+function load_controller!(integrator, ::PIController, file)
+    if !("time_integrator_qold" in keys(attributes(file)))
+        error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+    end
+    integrator.qold = read(attributes(file)["time_integrator_qold"])
+end
+
+# OrdinaryDiffEqCore < v4, PID controller:
+# Accumulated dt factor stored as integrator.qold.
+# Error history field layout changed between v3.31 and v3.32:
+#   <= v3.31: err::MVector{3} (single array field)
+#   v3.32–v3.33: err1, err2, err3 (individual scalar fields)
+function load_controller!(integrator, controller::PIDController, file)
+    if !("time_integrator_qold" in keys(attributes(file)) ||
+         !("time_integrator_controller_err" in keys(attributes(file))))
+        error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+    end
+    integrator.qold = read(attributes(file)["time_integrator_qold"])
+    if hasproperty(controller, :err)
+        # OrdinaryDiffEqCore <= v3.31
+        controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
+    elseif hasproperty(controller, :err1) && hasproperty(controller, :err2) &&
+           hasproperty(controller, :err3)
+        # OrdinaryDiffEqCore v3.32–v3.33
+        err = read(attributes(file)["time_integrator_controller_err"])
+        controller.err1 = err[1]
+        controller.err2 = err[2]
+        controller.err3 = err[3]
+    end
+end
+
+@static if Base.pkgversion(OrdinaryDiffEqCore) >= v"4"
+    # OrdinaryDiffEqCore >= v4, PI controller:
+    # Previous error estimate stored as PIControllerCache.errold (integrator.qold removed).
+    function load_controller!(integrator, controller::PIControllerCache, file)
+        if !("time_integrator_qold" in keys(attributes(file)))
+            error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+        end
+        controller.errold = read(attributes(file)["time_integrator_qold"])
+    end
+
+    # OrdinaryDiffEqCore >= v4, PID controller:
+    # Accumulated dt factor stored as PIDControllerCache.dt_factor (integrator.qold removed).
+    # Error history stored as PIDControllerCache.err::Vector{3}.
+    function load_controller!(integrator, controller::PIDControllerCache, file)
+        if !("time_integrator_qold" in keys(attributes(file)) ||
+             !("time_integrator_controller_err" in keys(attributes(file))))
+            error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+        end
+        controller.dt_factor = read(attributes(file)["time_integrator_qold"])
+        controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
+    end
+end
+
+# OrdinaryDiffEqCore < v4, PI controller:
+# Previous error estimate stored as integrator.qold.
+function store_controller!(file, ::PIController, integrator)
+    attributes(file)["time_integrator_qold"] = integrator.qold
+end
+
+# OrdinaryDiffEqCore < v4, PID controller:
+# Accumulated dt factor stored as integrator.qold.
+# Error history field layout changed between v3.31 and v3.32:
+#   <= v3.31: err::MVector{3} (single array field)
+#   v3.32–v3.33: err1, err2, err3 (individual scalar fields)
+function store_controller!(file, controller::PIDController, integrator)
+    attributes(file)["time_integrator_qold"] = integrator.qold
+    if hasproperty(controller, :err)
+        # OrdinaryDiffEqCore <= v3.31
+        attributes(file)["time_integrator_controller_err"] = controller.err
+    else
+        # OrdinaryDiffEqCore v3.32–v3.33
+        attributes(file)["time_integrator_controller_err"] = [controller.err1,
+            controller.err2,
+            controller.err3]
+    end
+end
+
+@static if Base.pkgversion(OrdinaryDiffEqCore) >= v"4"
+    # OrdinaryDiffEqCore >= v4, PI controller:
+    # Previous error estimate stored as PIControllerCache.errold (integrator.qold removed).
+    function store_controller!(file, controller::PIControllerCache, integrator)
+        attributes(file)["time_integrator_qold"] = controller.errold
+    end
+
+    # OrdinaryDiffEqCore >= v4, PID controller:
+    # Accumulated dt factor stored as PIDControllerCache.dt_factor (integrator.qold removed).
+    # Error history stored as PIDControllerCache.err::Vector{3}.
+    function store_controller!(file, controller::PIDControllerCache, integrator)
+        attributes(file)["time_integrator_qold"] = controller.dt_factor
+        attributes(file)["time_integrator_controller_err"] = controller.err
+    end
+end
+
+end

src/callbacks_step/callbacks_step.jl

@@ -30,6 +30,20 @@ end
            integrator.iter == integrator.opts.maxiters
 end
 
+# OrdinaryDiffEqCore < v4: legacy controller lives in integrator.opts.controller and
+# doubles as its own cache (setup_controller_cache returns the controller itself).
+# OrdinaryDiffEqCore >= v4: controller cache lives in integrator.controller_cache.
+get_controller_cache(integrator) = hasproperty(integrator, :controller_cache) ?
+                                   integrator.controller_cache :
+                                   integrator.opts.controller
+
+# OrdinaryDiffEq v7+ wraps the pure controller (parameters only) in controller_cache.controller;
+# older versions store it directly in controller_cache or opts.controller.
+function get_controller(integrator)
+    cc = get_controller_cache(integrator)
+    return hasproperty(cc, :controller) ? cc.controller : cc
+end
+
 # `include` callback definitions in the order that we currently prefer
 # when combining them into a `CallbackSet` which is called *after* a complete step
 # The motivation is as follows: The first callbacks belong to the current time step iteration:

src/callbacks_step/save_restart.jl

@@ -96,8 +96,7 @@ function (restart_callback::SaveRestartCallback)(integrator)
         save_restart_file(u_ode, t, dt, iter, semi, restart_callback)
         # If using an adaptive time stepping scheme, store controller values for restart
         if integrator.opts.adaptive
-            save_adaptive_time_integrator(integrator, integrator.opts.controller,
-                                          restart_callback)
+            save_adaptive_time_integrator(integrator, restart_callback)
         end
     end
 
@@ -170,29 +169,37 @@ Load the context information for time integrators with error-based step size con
 saved in a `restart_file`.
 """
 function load_adaptive_time_integrator!(integrator, restart_file::AbstractString)
-    controller = integrator.opts.controller
     # Read context information for controller
     h5open(restart_file, "r") do file
         # Ensure that the necessary information was saved
-        if !("time_integrator_qold" in keys(attributes(file))) ||
-           !("time_integrator_dtpropose" in keys(attributes(file))) ||
-           (hasproperty(controller, :err) &&
-            !("time_integrator_controller_err" in keys(attributes(file))))
+        if !("time_integrator_dtpropose" in keys(attributes(file)))
             error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
         end
-        # Load data that is required both for PIController and PIDController
-        integrator.qold = read(attributes(file)["time_integrator_qold"])
         integrator.dtpropose = read(attributes(file)["time_integrator_dtpropose"])
         # Accept step to use dtpropose already in the first step
         integrator.accept_step = true
         # Reevaluate integrator.fsal_first on the first step
         integrator.reeval_fsal = true
-        # Load additional parameters for PIDController
-        if hasproperty(controller, :err) # Distinguish PIDController from PIController
-            controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
-        end
+
+        load_controller!(integrator, file)
     end
 end
 
+function load_controller!(integrator, file)
+    return load_controller!(integrator, get_controller_cache(integrator), file)
+end
+
+function load_controller!(integrator, controller, file)
+    return error("Loading of controller $(typeof(controller)) not implemented.")
+end
+
+function store_controller!(file, integrator)
+    return store_controller!(file, get_controller_cache(integrator), integrator)
+end
+
+function store_controller!(file, controller, integrator)
+    return error("Storing of controller $(typeof(controller)) not implemented.")
+end
+
 include("save_restart_dg.jl")
 end # @muladd

src/callbacks_step/save_restart_dg.jl

@@ -336,8 +336,7 @@ function load_restart_file_on_root(mesh::Union{TreeMeshParallel, P4estMeshParall
 end
 
 # Store controller values for an adaptive time stepping scheme
-function save_adaptive_time_integrator(integrator,
-                                       controller, restart_callback)
+function save_adaptive_time_integrator(integrator, restart_callback)
     # Save only on root
     if mpi_isroot()
         @unpack output_directory = restart_callback
@@ -348,13 +347,8 @@ function save_adaptive_time_integrator(integrator,
 
         # Open file (preserve existing content)
         h5open(filename, "r+") do file
-            # Add context information as attributes both for PIController and PIDController
-            attributes(file)["time_integrator_qold"] = integrator.qold
             attributes(file)["time_integrator_dtpropose"] = integrator.dtpropose
-            # For PIDController is necessary to save additional parameters
-            if hasproperty(controller, :err) # Distinguish PIDController from PIController
-                attributes(file)["time_integrator_controller_err"] = controller.err
-            end
+            store_controller!(file, integrator)
         end
     end
 end

src/callbacks_step/save_solution_dg.jl

@@ -113,7 +113,7 @@ function save_solution_file(u, time, dt, timestep,
         data = u
         n_vars = nvariables(equations)
     else
-        data = map(u_node -> solution_variables(u_node, equations), u)
+        data = map(u_node -> solution_variables(u_node, equations), parent(u))
         # Find out variable count by looking at output from `solution_variables` function.
         n_vars = length(data[1])
     end

src/callbacks_step/summary.jl

@@ -208,7 +208,7 @@ function initialize_summary_callback(cb::DiscreteCallback, u, t, integrator;
         push!(setup,
               "abstol" => integrator.opts.abstol,
               "reltol" => integrator.opts.reltol,
-              "controller" => integrator.opts.controller)
+              "controller" => get_controller(integrator))
     end
     summary_box(io, "Time integration", setup)
     println()

src/postprocessing/spectral_analysis_2d.jl

@@ -9,7 +9,7 @@
     compute_kinetic_energy_spectrum(v1, v2)
 
 Compute an isotropic 1D kinetic energy spectrum from two 2D Cartesian velocity
-components `v1` and `v2`. For compressible Euler kinetic energy spectra, 
+components `v1` and `v2`. For compressible Euler kinetic energy spectra,
 pass density-weighted components `sqrt(rho) * v1` and `sqrt(rho) * v2`.
 The modal energy is normalized by `1 / N^2'.
 """
@@ -98,28 +98,28 @@ function interpolate_lgl_to_uniform_cartesian(u, mesh::TreeMesh{2},
     #      ↑     |             |
     #      |     |             |
     #      |     |             |
-    #      |     |             |   
+    #      |     |             |
     #      j    (1,1)--------(2,1)
     #             ---------> i (ξ)
     #
     #   Each of the 4 tree cells pastes a 2×2 block into the 4×4 array global grid, 'u_uniform'.
     #   Labels (column, row) are indices into that array: column ↔ x, row ↔ y on the domain.
     #   The assignment `u_uniform[:, r1, r2] .= interpolated` fills one such rectangle.
-    #   The bottom-right element of the original 2×2 block is labeled as 'E' in the diagram below and fills 
+    #   The bottom-right element of the original 2×2 block is labeled as 'E' in the diagram below and fills
     #   rows 3 to 4 and columns 3 to 4 of the global grid.
     #
     #         column  1         2         3         4
     #                 +---------+---------+---------+---------+
-    #         row 1   | (1,1)   | (2,1)   | (3,1)   | (4,1)   |  
-    #                 +---------+---------+---------+---------+   
-    #         row 2   | (1,2)   | (2,2)   | (3,2)   | (4,2)   |   
+    #         row 1   | (1,1)   | (2,1)   | (3,1)   | (4,1)   |
+    #                 +---------+---------+---------+---------+
+    #         row 2   | (1,2)   | (2,2)   | (3,2)   | (4,2)   |
     #                 +---------+---------+---------+---------+   ↑
-    #         row 3   | (1,3)   | (2,3)   | (3,3) E | (4,3) E |   |  
+    #         row 3   | (1,3)   | (2,3)   | (3,3) E | (4,3) E |   |
     #                 +---------+---------+---------+---------+   |  r2 = 3:4
     #         row 4   | (1,4)   | (2,4)   | (3,4) E | (4,4) E |   |
     #                 +---------+---------+---------+---------+   ↓
     #                                     |<-----r1 = 3:4---->|
-    #                 
+    #
 
     for element in eachelement(solver, cache)
         # Gather conservative nodal values on the reference LGL tensor grid for the element
@@ -168,7 +168,7 @@ function compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{2},
                                          equations::AbstractCompressibleEulerEquations,
                                          dg::DGMultiSBP, cache)
     # Unpacks the primitive variables from the conservative state for FDSBP DGMulti solutions
-    u_values = StructArray(u)
+    u_values = parent(u)
     n_points = length(u_values)
     n = round(Int, sqrt(n_points))
     q = cons2prim.(u_values, Ref(equations)) # q is the vector that contains the primitive variables for density and velocity converted from the conservative variables

src/postprocessing/spectral_analysis_3d.jl

@@ -9,7 +9,7 @@
     compute_kinetic_energy_spectrum(v1, v2, v3)
 
 Compute an isotropic 1D kinetic energy spectrum from three 3D Cartesian velocity
-components `v1`, `v2`, `v3`. For compressible Euler kinetic energy spectra, 
+components `v1`, `v2`, `v3`. For compressible Euler kinetic energy spectra,
 pass density-weighted components `sqrt(rho) * v1`, `sqrt(rho) * v2`, and `sqrt(rho) * v3`.
 The modal energy is normalized by `1 / N^3`.
 """
@@ -145,7 +145,7 @@ function compute_kinetic_energy_spectrum(u, mesh::DGMultiMesh{3},
                                          equations::AbstractCompressibleEulerEquations,
                                          dg::DGMultiSBP, cache)
     # Unpacks the primitive variables from the conservative state for FDSBP DGMulti solutions
-    u_values = StructArray(u)
+    u_values = parent(u)
     n_points = length(u_values)
     n = round(Int, n_points^(1 / 3))
     q = cons2prim.(u_values, Ref(equations)) # q is the vector that contains the primitive variables for density and velocity converted from the conservative variables