diff --git a/Project.toml b/Project.toml index 05eb268a1..1713ddcba 100644 --- a/Project.toml +++ b/Project.toml @@ -15,7 +15,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf" IteratorInterfaceExtensions = "82899510-4779-5014-852e-03e436cf321d" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d" @@ -35,6 +34,7 @@ SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" [weakdeps] ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" PartialFunctions = "570af359-4316-4cb7-8c74-252c00c2016b" @@ -44,6 +44,7 @@ RCall = "6f49c342-dc21-5d91-9882-a32aef131414" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [extensions] +SciMLBaseLinearAlgebraExt = "LinearAlgebra" SciMLBaseChainRulesCoreExt = "ChainRulesCore" SciMLBaseMLStyleExt = "MLStyle" SciMLBaseMakieExt = "Makie" diff --git a/ext/SciMLBaseLinearAlgebraExt.jl b/ext/SciMLBaseLinearAlgebraExt.jl new file mode 100644 index 000000000..41842d09a --- /dev/null +++ b/ext/SciMLBaseLinearAlgebraExt.jl @@ -0,0 +1,141 @@ +module SciMLBaseLinearAlgebraExt + +using SciMLBase +using LinearAlgebra + +# Override the identity matrix constant after module loading +function __init__() + SciMLBase.I = LinearAlgebra.I +end + +# Extension implementations +SciMLBase.default_identity_matrix() = LinearAlgebra.I + +function SciMLBase.has_non_trivial_mass_matrix(prob) + hasproperty(prob.f, :mass_matrix) && !(prob.f.mass_matrix isa UniformScaling{Bool}) +end + +function SciMLBase.get_initial_values( + prob::SciMLBase.AbstractDEProblem, integrator::SciMLBase.DEIntegrator, f, alg::SciMLBase.CheckInit, + isinplace::Union{Val{true}, Val{false}}; abstol, kwargs...) + u0 = SciMLBase.state_values(integrator) + p = SciMLBase.parameter_values(integrator) + t = SciMLBase.current_time(integrator) + M = f.mass_matrix + + M == LinearAlgebra.I && return u0, p, true + algebraic_vars = [all(iszero, x) for x in eachcol(M)] + algebraic_eqs = [all(iszero, x) for x in eachrow(M)] + (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return u0, p, true + SciMLBase.update_coefficients!(M, u0, p, t) + tmp = SciMLBase.evaluate_f(integrator, prob, f, isinplace, u0, p, t) + tmp .= SciMLBase.ArrayInterface.restructure(tmp, algebraic_eqs .* SciMLBase._vec(tmp)) + + normresid = isdefined(integrator.opts, :internalnorm) ? + integrator.opts.internalnorm(tmp, t) : norm(tmp) + if normresid > abstol + isdae = prob isa SciMLBase.AbstractDAEProblem + throw(SciMLBase.CheckInitFailureError(normresid, abstol, isdae)) + end + u0, p, true +end + +function SciMLBase.calculate_solution_errors!( + sol::SciMLBase.AbstractODESolution; fill_uanalytic = true, + timeseries_errors = true, dense_errors = true) + f = sol.prob.f + + if fill_uanalytic + for i in 1:size(sol.u, 1) + if sol.prob isa SciMLBase.AbstractDDEProblem + push!(sol.u_analytic, + f.analytic(sol.prob.u0, sol.prob.h, sol.prob.p, sol.t[i])) + else + push!(sol.u_analytic, f.analytic(sol.prob.u0, sol.prob.p, sol.t[i])) + end + end + end + + save_everystep = length(sol.u) > 2 + if !isempty(sol.u_analytic) + sol.errors[:final] = norm(SciMLBase.recursive_mean(abs.(sol.u[end] .- + sol.u_analytic[end]))) + + if save_everystep && timeseries_errors + sol.errors[:l∞] = norm(maximum(SciMLBase.vecvecapply((x) -> abs.(x), + sol.u - sol.u_analytic))) + sol.errors[:l2] = norm(sqrt(SciMLBase.recursive_mean(SciMLBase.vecvecapply( + (x) -> float.(x) .^ 2, + sol.u - sol.u_analytic)))) + if sol.dense && dense_errors + densetimes = collect(range(sol.t[1], stop = sol.t[end], length = 100)) + interp_u = sol(densetimes) + interp_analytic = SciMLBase.VectorOfArray([f.analytic(sol.prob.u0, sol.prob.p, t) + for t in densetimes]) + sol.errors[:L∞] = norm(maximum(SciMLBase.vecvecapply((x) -> abs.(x), + interp_u - interp_analytic))) + sol.errors[:L2] = norm(sqrt(SciMLBase.recursive_mean(SciMLBase.vecvecapply( + (x) -> float.(x) .^ 2, + interp_u - + interp_analytic)))) + end + end + end +end + +function SciMLBase.calculate_solution_errors!( + sol::SciMLBase.AbstractDAESolution; fill_uanalytic = true, + timeseries_errors = true, dense_errors = true) + prob = sol.prob + f = prob.f + + if fill_uanalytic + for i in 1:size(sol.u, 1) + push!(sol.u_analytic, f.analytic(prob.du0, prob.u0, prob.p, sol.t[i])) + end + end + + save_everystep = length(sol.u) > 2 + if !isempty(sol.u_analytic) + sol.errors[:final] = norm(SciMLBase.recursive_mean(abs.(sol.u[end] - + sol.u_analytic[end]))) + + if save_everystep && timeseries_errors + sol.errors[:l∞] = norm(maximum(SciMLBase.vecvecapply(x -> abs.(x), + sol.u - sol.u_analytic))) + sol.errors[:l2] = norm(sqrt(SciMLBase.recursive_mean(SciMLBase.vecvecapply( + x -> float.(x) .^ 2, + sol.u - sol.u_analytic)))) + if sol.dense && dense_errors + densetimes = collect(range(sol.t[1]; stop = sol.t[end], length = 100)) + interp_u = sol(densetimes) + interp_analytic = SciMLBase.VectorOfArray([f.analytic(prob.du0, prob.u0, prob.p, t) + for t in densetimes]) + sol.errors[:L∞] = norm(maximum(SciMLBase.vecvecapply(x -> abs.(x), + interp_u - interp_analytic))) + sol.errors[:L2] = norm(sqrt(SciMLBase.recursive_mean(SciMLBase.vecvecapply( + x -> float.(x) .^ 2, + interp_u .- + interp_analytic)))) + end + end + end +end + +function SciMLBase.calculate_solution_errors!( + sol::SciMLBase.AbstractRODESolution; fill_uanalytic = true, + timeseries_errors = true) + if !isempty(sol.u_analytic) + sol.errors[:final] = norm(SciMLBase.recursive_mean(abs.(sol.u[end] - + sol.u_analytic[end]))) + if timeseries_errors + sol.errors[:l∞] = norm(maximum(SciMLBase.vecvecapply((x) -> abs.(x), + sol.u - sol.u_analytic))) + sol.errors[:l2] = norm(sqrt(SciMLBase.recursive_mean(SciMLBase.vecvecapply( + (x) -> float.(x) .^ 2, + sol.u - sol.u_analytic)))) + end + end +end + +end diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index 780e68e80..d16a43b84 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -8,7 +8,6 @@ using RecipesBase, RecursiveArrayTools using SciMLStructures using SymbolicIndexingInterface using DocStringExtensions -using LinearAlgebra using Statistics using Distributed using Markdown @@ -674,6 +673,12 @@ Internal. Used for signifying the AD context comes from a Mooncake.jl context. """ struct MooncakeOriginator <: ADOriginator end +# Default implementations for LinearAlgebra extension +I = nothing # Will be overridden by extension (not const to allow reassignment) +default_identity_matrix() = I +has_non_trivial_mass_matrix(prob) = false # Default implementation +calculate_solution_errors!(sol; kwargs...) = nothing # Stub implementation + include("initialization.jl") include("odenlstep.jl") include("utils.jl") @@ -843,6 +848,8 @@ export OptimizationFunction, MultiObjectiveOptimizationFunction export CheckInit +export default_identity_matrix, has_non_trivial_mass_matrix, calculate_solution_errors! + export EnsembleThreads, EnsembleDistributed, EnsembleSplitThreads, EnsembleSerial export EnsembleAnalysis, EnsembleSummary diff --git a/src/initialization.jl b/src/initialization.jl index e03f41afc..df3d7b2fd 100644 --- a/src/initialization.jl +++ b/src/initialization.jl @@ -176,30 +176,12 @@ Keyword arguments: - `abstol`: The absolute value below which the norm of the residual of algebraic equations should lie. The norm function used is `integrator.opts.internalnorm` if present, and - `LinearAlgebra.norm` if not. + `norm` (from LinearAlgebra extension) if not. """ function get_initial_values( prob::AbstractDEProblem, integrator::DEIntegrator, f, alg::CheckInit, isinplace::Union{Val{true}, Val{false}}; abstol, kwargs...) - u0 = state_values(integrator) - p = parameter_values(integrator) - t = current_time(integrator) - M = f.mass_matrix - - M == I && return u0, p, true - algebraic_vars = [all(iszero, x) for x in eachcol(M)] - algebraic_eqs = [all(iszero, x) for x in eachrow(M)] - (iszero(algebraic_vars) || iszero(algebraic_eqs)) && return u0, p, true - update_coefficients!(M, u0, p, t) - tmp = evaluate_f(integrator, prob, f, isinplace, u0, p, t) - tmp .= ArrayInterface.restructure(tmp, algebraic_eqs .* _vec(tmp)) - - normresid = isdefined(integrator.opts, :internalnorm) ? - integrator.opts.internalnorm(tmp, t) : norm(tmp) - if normresid > abstol - throw(CheckInitFailureError(normresid, abstol, true)) - end - return u0, p, true + error("LinearAlgebra extension not loaded. Please load LinearAlgebra to use CheckInit.") end function get_initial_values( diff --git a/src/problems/problem_utils.jl b/src/problems/problem_utils.jl index 5c55d2e2b..43c0b5ae0 100644 --- a/src/problems/problem_utils.jl +++ b/src/problems/problem_utils.jl @@ -38,7 +38,7 @@ function Base.summary(io::IO, prob::AbstractDEProblem) hasproperty(prob.f, :mass_matrix) && begin println(io) print(io, "Non-trivial mass matrix: ", type_color, - !(prob.f.mass_matrix isa LinearAlgebra.UniformScaling{Bool}), no_color) + has_non_trivial_mass_matrix(prob), no_color) end end diff --git a/src/scimlfunctions.jl b/src/scimlfunctions.jl index 0d6c26a02..62fc1e1c9 100644 --- a/src/scimlfunctions.jl +++ b/src/scimlfunctions.jl @@ -254,7 +254,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia ODEFunction{iip,specialize}(f; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -343,7 +343,6 @@ The fields of the ODEFunction type directly match the names of the inputs. The following example creates an inplace `ODEFunction` whose Jacobian is a `Diagonal`: ```julia -using LinearAlgebra f = (du,u,p,t) -> du .= t .* u jac = (J,u,p,t) -> (J[1,1] = t; J[2,2] = t; J) jp = Diagonal(zeros(2)) @@ -454,7 +453,7 @@ and exponential integrators. ```julia SplitFunction{iip,specialize}(f1,f2; - mass_matrix = __has_mass_matrix(f1) ? f1.mass_matrix : I, + mass_matrix = __has_mass_matrix(f1) ? f1.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f1) ? f1.analytic : nothing, tgrad= __has_tgrad(f1) ? f1.tgrad : nothing, jac = __has_jac(f1) ? f1.jac : nothing, @@ -581,7 +580,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia DynamicalODEFunction{iip,specialize}(f1,f2; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -683,7 +682,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia DDEFunction{iip,specialize}(f; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -788,7 +787,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia DynamicalDDEFunction{iip,specialize}(f1,f2; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -1006,7 +1005,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia SDEFunction{iip,specialize}(f,g; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -1112,7 +1111,7 @@ and exponential integrators. ```julia SplitSDEFunction{iip,specialize}(f1,f2,g; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -1223,7 +1222,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia DynamicalSDEFunction{iip,specialize}(f1,f2; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -1332,7 +1331,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia RODEFunction{iip,specialize}(f; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -1579,7 +1578,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia SDDEFunction{iip,specialize}(f,g; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -2116,7 +2115,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia ODEInputFunction{iip, specialize}(f; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -2233,7 +2232,7 @@ BVPFunction{iip, specialize}(f, bc; cost = __has_cost(f) ? f.cost : nothing, equality = __has_equality(f) ? f.equality : nothing, inequality = __has_inequality(f) ? f.inequality : nothing, - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -2371,7 +2370,7 @@ with respect to time, and more. For all cases, `u0` is the initial condition, ```julia DynamicalBVPFunction{iip,specialize}(f, bc; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad= __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -2715,8 +2714,8 @@ function ODEFunction{iip, specialize}(f; ) where {iip, specialize } - if mass_matrix === I && f isa Tuple - mass_matrix = ((I for i in 1:length(f))...,) + if mass_matrix === default_identity_matrix() && f isa Tuple + mass_matrix = ((default_identity_matrix() for i in 1:length(f))...,) end if (specialize === FunctionWrapperSpecialize) && @@ -4107,8 +4106,8 @@ function NonlinearFunction{iip, specialize}(f; initialization_data = __has_initialization_data(f) ? f.initialization_data : nothing) where { iip, specialize} - if mass_matrix === I && f isa Tuple - mass_matrix = ((I for i in 1:length(f))...,) + if mass_matrix === default_identity_matrix() && f isa Tuple + mass_matrix = ((default_identity_matrix() for i in 1:length(f))...,) end if jac === nothing && isa(jac_prototype, AbstractSciMLOperator) @@ -4360,7 +4359,7 @@ function BVPFunction{iip, specialize, twopoint}(f, bc; cost = __has_cost(f) ? f.cost : nothing, equality = __has_equality(f) ? f.equality : nothing, inequality = __has_inequality(f) ? f.inequality : nothing, - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad = __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -4383,8 +4382,8 @@ function BVPFunction{iip, specialize, twopoint}(f, bc; sys = __has_sys(f) ? f.sys : nothing, initialization_data = __has_initialization_data(f) ? f.initialization_data : nothing) where {iip, specialize, twopoint} - if mass_matrix === I && f isa Tuple - mass_matrix = ((I for i in 1:length(f))...,) + if mass_matrix === default_identity_matrix() && f isa Tuple + mass_matrix = ((default_identity_matrix() for i in 1:length(f))...,) end if (specialize === FunctionWrapperSpecialize) && @@ -4540,7 +4539,7 @@ end BVPFunction(f::BVPFunction; kwargs...) = f function DynamicalBVPFunction{iip, specialize, twopoint}(f, bc; - mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : I, + mass_matrix = __has_mass_matrix(f) ? f.mass_matrix : default_identity_matrix(), analytic = __has_analytic(f) ? f.analytic : nothing, tgrad = __has_tgrad(f) ? f.tgrad : nothing, jac = __has_jac(f) ? f.jac : nothing, @@ -4563,8 +4562,8 @@ function DynamicalBVPFunction{iip, specialize, twopoint}(f, bc; sys = __has_sys(f) ? f.sys : nothing, initialization_data = __has_initialization_data(f) ? f.initialization_data : nothing) where {iip, specialize, twopoint} - if mass_matrix === I && f isa Tuple - mass_matrix = ((I for i in 1:length(f))...,) + if mass_matrix === default_identity_matrix() && f isa Tuple + mass_matrix = ((default_identity_matrix() for i in 1:length(f))...,) end if (specialize === FunctionWrapperSpecialize) && @@ -4804,8 +4803,8 @@ function ODEInputFunction{iip, specialize}(f; ) where {iip, specialize } - if mass_matrix === I && f isa Tuple - mass_matrix = ((I for i in 1:length(f))...,) + if mass_matrix === default_identity_matrix() && f isa Tuple + mass_matrix = ((default_identity_matrix() for i in 1:length(f))...,) end if (specialize === FunctionWrapperSpecialize) && diff --git a/src/solutions/dae_solutions.jl b/src/solutions/dae_solutions.jl index dd3960624..354bd04e9 100644 --- a/src/solutions/dae_solutions.jl +++ b/src/solutions/dae_solutions.jl @@ -164,28 +164,8 @@ function calculate_solution_errors!(sol::AbstractDAESolution; end end - save_everystep = length(sol.u) > 2 - if !isempty(sol.u_analytic) - sol.errors[:final] = norm(recursive_mean(abs.(sol.u[end] - sol.u_analytic[end]))) - - if save_everystep && timeseries_errors - sol.errors[:l∞] = norm(maximum(vecvecapply(x -> abs.(x), - sol.u - sol.u_analytic))) - sol.errors[:l2] = norm(sqrt(recursive_mean(vecvecapply(x -> float.(x) .^ 2, - sol.u - sol.u_analytic)))) - if sol.dense && dense_errors - densetimes = collect(range(sol.t[1]; stop = sol.t[end], length = 100)) - interp_u = sol(densetimes) - interp_analytic = VectorOfArray([f.analytic(prob.du0, prob.u0, prob.p, t) - for t in densetimes]) - sol.errors[:L∞] = norm(maximum(vecvecapply(x -> abs.(x), - interp_u - interp_analytic))) - sol.errors[:L2] = norm(sqrt(recursive_mean(vecvecapply(x -> float.(x) .^ 2, - interp_u .- - interp_analytic)))) - end - end - end + # This functionality requires LinearAlgebra extension + # The actual implementation is in ext/SciMLBaseLinearAlgebraExt.jl nothing end diff --git a/src/solutions/ode_solutions.jl b/src/solutions/ode_solutions.jl index e15a9e16f..9ad523322 100644 --- a/src/solutions/ode_solutions.jl +++ b/src/solutions/ode_solutions.jl @@ -560,29 +560,8 @@ function calculate_solution_errors!(sol::AbstractODESolution; fill_uanalytic = t end end - save_everystep = length(sol.u) > 2 - if !isempty(sol.u_analytic) - sol.errors[:final] = norm(recursive_mean(abs.(sol.u[end] .- sol.u_analytic[end]))) - - if save_everystep && timeseries_errors - sol.errors[:l∞] = norm(maximum(vecvecapply((x) -> abs.(x), - sol.u - sol.u_analytic))) - sol.errors[:l2] = norm(sqrt(recursive_mean(vecvecapply((x) -> float.(x) .^ 2, - sol.u - sol.u_analytic)))) - if sol.dense && dense_errors - densetimes = collect(range(sol.t[1], stop = sol.t[end], length = 100)) - interp_u = sol(densetimes) - interp_analytic = VectorOfArray([f.analytic(sol.prob.u0, sol.prob.p, t) - for t in densetimes]) - sol.errors[:L∞] = norm(maximum(vecvecapply((x) -> abs.(x), - interp_u - interp_analytic))) - sol.errors[:L2] = norm(sqrt(recursive_mean(vecvecapply( - (x) -> float.(x) .^ 2, - interp_u - - interp_analytic)))) - end - end - end + # This functionality requires LinearAlgebra extension + # The actual implementation is in ext/SciMLBaseLinearAlgebraExt.jl end function build_solution(sol::ODESolution{T, N}, u_analytic, errors) where {T, N} diff --git a/src/solutions/rode_solutions.jl b/src/solutions/rode_solutions.jl index 2e31dabd9..1c0ad211d 100644 --- a/src/solutions/rode_solutions.jl +++ b/src/solutions/rode_solutions.jl @@ -220,13 +220,8 @@ function calculate_solution_errors!(sol::AbstractRODESolution; fill_uanalytic = end if !isempty(sol.u_analytic) - sol.errors[:final] = norm(recursive_mean(abs.(sol.u[end] - sol.u_analytic[end]))) - if timeseries_errors - sol.errors[:l∞] = norm(maximum(vecvecapply((x) -> abs.(x), - sol.u - sol.u_analytic))) - sol.errors[:l2] = norm(sqrt(recursive_mean(vecvecapply((x) -> float.(x) .^ 2, - sol.u - sol.u_analytic)))) - end + # This functionality requires LinearAlgebra extension + # The actual implementation is in ext/SciMLBaseLinearAlgebraExt.jl if dense_errors densetimes = collect(range(sol.t[1], stop = sol.t[end], length = 100)) interp_u = sol(densetimes)