Fix ChunkSplitters v3 API (precompile failure), re-export solve, add stdlib compat

Project.toml

@@ -1,7 +1,7 @@
 name = "FiniteVolumeMethod"
 uuid = "d4f04ab7-4f65-4d72-8a28-7087bc7f46f4"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.2.1"
+version = "1.2.2"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
@@ -14,12 +14,14 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-ChunkSplitters = "0.1, 1.0, 2.0, 3.2"
-CommonSolve = "0.2"
-DelaunayTriangulation = "1.0"
-PreallocationTools = "0.4, 1.2"
-PrecompileTools = "1.2"
-SciMLBase = "2.34, 3.1"
+ChunkSplitters = "3.2"
+CommonSolve = "0.2.4"
+DelaunayTriangulation = "1.6.6"
+LinearAlgebra = "1"
+PreallocationTools = "1.2"
+PrecompileTools = "1.2.1"
+SciMLBase = "3.1"
+SparseArrays = "1"
 Test = "1"
 julia = "1"
 

README.md

@@ -42,7 +42,7 @@ sol = solve(prob, Tsit5(), saveat = 0.001)
 u = Observable(sol.u[1])
 fig, ax, sc = tricontourf(tri, u, levels = 0:5:50, colormap = :matter)
 tightlimits!(ax)
-record(fig, "anim.gif", eachindex(sol)) do i
+record(fig, "anim.gif", eachindex(sol.u)) do i
     u[] = sol.u[i]
 end
 ```

docs/Project.toml

@@ -25,11 +25,13 @@ LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 PProf = "e4faabce-9ead-11e9-39d9-4379958e3056"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SimpleGraphs = "55797a34-41de-5266-9ec1-32ac4eb504d3"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
@@ -64,11 +66,13 @@ LiveServer = "1.5.0"
 NaturalNeighbours = "1.3.6"
 NonlinearSolve = "4.19.1"
 OrdinaryDiffEq = "7.0.0"
+OrdinaryDiffEqSDIRK = "2.7.0"
 PProf = "3.2.0"
 ReferenceTests = "0.10.6"
 Roots = "3.0.0"
 SciMLBase = "3.18.0"
-SimpleGraphs = "0.8.10"
+SciMLOperators = "1.22.0"
+SimpleGraphs = "0.8.8 - 0.8"
 Sparspak = "0.3.15"
 SpecialFunctions = "2.8.0"
 StableRNGs = "1.0.4"

docs/src/literate_tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.jl

@@ -79,13 +79,13 @@ flux = (x, y, t, α, β, γ, p) -> (-α, -β)
 # We now solve the problem. We provide the solver for this problem.
 # In my experience, I've found that `TRBDF2(linsolve=KLUFactorization())` typically
 # has the best performance for these problems.
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.01, parallel = Val(false))
 ind = findall(DelaunayTriangulation.each_point_index(tri)) do i #hide
     !DelaunayTriangulation.has_vertex(tri, i) #hide
 end #hide
 using Test #hide
-@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol)), :, length(sol))[ind, :] # make sure that missing vertices don't change #hide
+@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol.u)), :, length(sol.u))[ind, :] # make sure that missing vertices don't change #hide
 sol |> tc #hide
 
 #-
@@ -151,12 +151,12 @@ function exact_solution(x, y, t, A, ζ, f, α) #src
 end #src
 function compare_solutions(sol, tri, α, f) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
     ζ = get_ζ_terms(20, 20, α) #src
     A = get_sum_coefficients(20, 20, α, ζ) #src
-    for i in eachindex(sol) #src
+    for i in eachindex(sol.u) #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
             u[j, i] = exact_solution(x[j, i], y[j, i], sol.t[i], A, ζ, f, α) #src
@@ -166,7 +166,7 @@ function compare_solutions(sol, tri, α, f) #src
 end #src
 x, y, u = compare_solutions(sol, tri, α, f) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     local ax #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 0:0.01:1, colormap = :matter) #src

docs/src/literate_tutorials/diffusion_equation_on_a_square_plate.jl

@@ -94,10 +94,10 @@ function exact_solution(x, y, t) #src
 end #src
 function compare_solutions(sol, tri) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
             u[j, i] = exact_solution(x[j, i], y[j, i], sol.t[i]) #src
@@ -107,7 +107,7 @@ function compare_solutions(sol, tri) #src
 end #src
 x, y, u = compare_solutions(sol, tri) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     local ax #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 0:5:50, colormap = :matter) #src

docs/src/literate_tutorials/diffusion_equation_on_an_annulus.jl

@@ -79,7 +79,7 @@ prob = FVMProblem(
 )
 
 #-
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.2)
 sol |> tc #hide
 

docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl

@@ -63,7 +63,7 @@ v_prob = FVMProblem(
 prob = FVMSystem(u_prob, v_prob)
 
 # Now that we have our system, we can solve.
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 10.0, parallel = Val(false))
 sol |> tc #hide
 
@@ -79,7 +79,7 @@ y = LinRange(-1, 1, 200)
 heatmap!(ax, x, y, u, colorrange = (0.0, 0.4))
 hidedecorations!(ax)
 record(
-    fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
+    fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol.u);
     framerate = 60
 ) do _i
     i[] = _i

docs/src/literate_tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.jl

@@ -64,7 +64,7 @@ steady_prob = SteadyFVMProblem(prob)
 using NonlinearSolve
 sol = solve(steady_prob, NewtonRaphson())
 copyto!(prob.initial_condition, sol.u) # this also changes steady_prob's initial condition
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 sol |> tc #hide
 

docs/src/literate_tutorials/laplaces_equation_with_internal_dirichlet_conditions.jl

@@ -116,7 +116,7 @@ prob = FVMProblem(
 steady_prob = SteadyFVMProblem(prob)
 
 # Now let's solve the problem.
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 sol |> tc #hide
 

docs/src/literate_tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.jl

@@ -25,7 +25,7 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # with \eqref{eq:advdiffeq}. For the mesh, we could use
 # `triangulate_rectangle`, but we want to put most of the triangles
 # near the origin, so we need to use `refine!` on an initial mesh.
-using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie
+using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie, StableRNGs
 L = 30
 tri = triangulate_rectangle(-L, L, -L, L, 2, 2, single_boundary = true)
 tot_area = get_area(tri)
@@ -39,7 +39,7 @@ area_constraint = (_tri, T) -> begin
     flag = A ≥ max_area_function(A, dist_to_origin)
     return flag
 end
-refine!(tri; min_angle = 33.0, custom_constraint = area_constraint)
+refine!(tri; min_angle = 33.0, custom_constraint = area_constraint, rng = StableRNG(123))
 triplot(tri)
 
 #-
@@ -96,7 +96,7 @@ prob = FVMProblem(
 )
 
 # Now we can solve and visualise the solution.
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 times = [0, 10, 25, 50, 100, 200, 250]
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = times)
 sol |> tc #hide
@@ -105,7 +105,7 @@ sol |> tc #hide
 using CairoMakie
 using ReferenceTests #src
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(
         fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
@@ -133,9 +133,9 @@ function exact_solution(x, y, t, D, ν) #src
     return 1 / (4 * D * π * t) * exp(1 / (4 * D * t) * (-(x - ν * t)^2 - y^2)) #src
 end #src
 function get_errs(_sol, tri, flux_parameters) #src
-    _errs = zeros(length(_sol)) #src
+    _errs = zeros(length(_sol.u)) #src
     _err = zeros(DelaunayTriangulation.num_points(tri)) #src
-    for i in eachindex(_sol) #src
+    for i in eachindex(_sol.u) #src
         !DelaunayTriangulation.has_vertex(tri, i) && continue #src
         i == 1 && continue #src
         m = maximum(_sol.u[i]) #src
@@ -205,8 +205,8 @@ for j in eachindex(y)
         triangles[i, j] = jump_and_march(tri, (x[i], y[j]))
     end
 end
-interpolated_vals = zeros(length(x), length(y), length(sol))
-for k in eachindex(sol)
+interpolated_vals = zeros(length(x), length(y), length(sol.u))
+for k in eachindex(sol.u)
     for j in eachindex(y)
         for i in eachindex(x)
             interpolated_vals[i, j, k] = pl_interpolate(
@@ -218,9 +218,9 @@ end
 
 # Let's visualise these results to check their accuracy. We compute the triangulation of
 # our grid to make the `tricontourf` call faster.
-_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]')
+_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]', rng = StableRNG(123))
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(
         fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",

docs/src/literate_tutorials/porous_fisher_equation_and_travelling_waves.jl

@@ -47,7 +47,7 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # takes the solution some time to evolve towards the travelling wave solution.
 #
 # Now with this preamble out of the way, let us solve this problem.
-using DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, LinearSolve
+using DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 a, b, c, d, nx, ny = 0.0, 3.0, 0.0, 40.0, 60, 80
 tri = triangulate_rectangle(a, b, c, d, nx, ny; single_boundary = false)
 mesh = FVMGeometry(tri)
@@ -81,10 +81,10 @@ c = sqrt(λ / (2D))
 cₘᵢₙ = sqrt(λ * D / 2)
 zᶜ = 0.0
 exact_solution(z) = ifelse(z ≤ zᶜ, 1 - exp(cₘᵢₙ * z), zero(z))
-travelling_wave_values = zeros(ny, length(sol) - large_time_idx + 1)
+travelling_wave_values = zeros(ny, length(sol.u) - large_time_idx + 1)
 z_vals = zero(travelling_wave_values)
 u_mat = [reshape(u, (nx, ny)) for u in sol.u]
-for (i, t_idx) in pairs(large_time_idx:lastindex(sol))
+for (i, t_idx) in pairs(large_time_idx:lastindex(sol.u))
     u = u_mat[t_idx]
     τ = sol.t[t_idx]
     for k in 1:ny
@@ -109,10 +109,10 @@ for (i, j) in zip(1:3, (1, 51, 101))
     tightlimits!(ax)
 end
 ax = Axis(fig[1, 4], width = 900, height = 600)
-colors = cgrad(:matter, length(sol) - large_time_idx + 1; categorical = false)
+colors = cgrad(:matter, length(sol.u) - large_time_idx + 1; categorical = false)
 [
     lines!(ax, z_vals[:, i], travelling_wave_values[:, i], color = colors[i], linewidth = 2)
-        for i in 1:(length(sol) - large_time_idx + 1)
+        for i in 1:(length(sol.u) - large_time_idx + 1)
 ]
 exact_z_vals = collect(LinRange(extrema(z_vals)..., 500))
 exact_travelling_wave_values = exact_solution.(exact_z_vals)

docs/src/literate_tutorials/porous_medium_equation.jl

@@ -60,7 +60,7 @@ prob = FVMProblem(
 
 #-
 ## Step 4: Solve
-using LinearSolve, OrdinaryDiffEq
+using LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()); saveat = 3.0)
 sol |> tc #hide
 
@@ -99,10 +99,10 @@ function exact_solution(x, y, t, m, M, D) #src
 end #src
 function compare_solutions(sol, tri, m, M, D) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         !DelaunayTriangulation.has_vertex(tri, i) && continue #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
@@ -113,7 +113,7 @@ function compare_solutions(sol, tri, m, M, D) #src
 end #src
 x, y, u = compare_solutions(sol, tri, m, M, D) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 0:0.005:0.1, colormap = :matter) #src
     ax = Axis(fig[2, i], width = 600, height = 600) #src
@@ -207,10 +207,10 @@ function exact_solution(x, y, t, m, M, D, λ) #src
 end #src
 function compare_solutions(sol, tri, m, M, D, λ) #src
     n = DelaunayTriangulation.num_solid_vertices(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
             u[j, i] = exact_solution(x[j, i], y[j, i], sol.t[i], m, M, D, λ) #src
@@ -220,7 +220,7 @@ function compare_solutions(sol, tri, m, M, D, λ) #src
 end #src
 x, y, u = compare_solutions(sol, tri, m, M, D, λ) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(
         ax, tri, sol.u[i], levels = 0:0.05:1, extendlow = :auto,

docs/src/literate_tutorials/reaction_diffusion_brusselator_system_of_pdes.jl

@@ -143,7 +143,7 @@ mesh = FVMGeometry(tri)
 system = FVMSystem(Φ_prob, Ψ_prob)
 
 # We can now solve the problem just as we've done previously.
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(system, TRBDF2(linsolve = KLUFactorization()), saveat = 1.0)
 sol |> tc #hide
 
@@ -159,7 +159,7 @@ sol.u[3][1, :] |> tc #hide
 # are the values of $\Psi$ at the third time. We can visualise the solutions as follows:
 using CairoMakie
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax1 = Axis(
         fig[1, i], xlabel = L"x", ylabel = L"y",
         width = 400, height = 400,
@@ -170,8 +170,8 @@ for i in eachindex(sol)
         width = 400, height = 400,
         title = L"\Psi: t = %$(sol.t[i])", titlealign = :left
     )
-    tricontourf!(ax1, tri, sol[i][1, :], levels = 0:0.1:1, colormap = :matter)
-    tricontourf!(ax2, tri, sol[i][2, :], levels = 1:10:100, colormap = :matter)
+    tricontourf!(ax1, tri, sol.u[i][1, :], levels = 0:0.1:1, colormap = :matter)
+    tricontourf!(ax2, tri, sol.u[i][2, :], levels = 1:10:100, colormap = :matter)
 end
 resize_to_layout!(fig)
 fig
@@ -180,7 +180,7 @@ using ReferenceTests #src
 
 x = getx.(get_points(tri)) #src
 y = gety.(get_points(tri)) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax3 = Axis(
         fig[3, i], xlabel = L"x", ylabel = L"y", #src
         width = 400, height = 400, #src