diff --git a/src/SciMLBase.jl b/src/SciMLBase.jl index 986ec0949..e930101e7 100644 --- a/src/SciMLBase.jl +++ b/src/SciMLBase.jl @@ -791,9 +791,10 @@ export isinplace export solve, solve!, init, discretize, symbolic_discretize -export LinearProblem, NonlinearProblem, IntervalNonlinearProblem, - IntegralProblem, SampledIntegralProblem, OptimizationProblem, - NonlinearLeastSquaresProblem +export LinearProblem, IntervalNonlinearProblem, + IntegralProblem, SampledIntegralProblem, OptimizationProblem + +export NonlinearProblem, SCCNonlinearProblem, NonlinearLeastSquaresProblem export DiscreteProblem, ImplicitDiscreteProblem export SteadyStateProblem, SteadyStateSolution diff --git a/src/problems/nonlinear_problems.jl b/src/problems/nonlinear_problems.jl index 78a6e05ef..58241f289 100644 --- a/src/problems/nonlinear_problems.jl +++ b/src/problems/nonlinear_problems.jl @@ -321,3 +321,148 @@ end function NonlinearLeastSquaresProblem(f, u0, p = NullParameters(); kwargs...) return NonlinearLeastSquaresProblem(NonlinearFunction(f), u0, p; kwargs...) end + +@doc doc""" + SCCNonlinearProblem(probs, explicitfuns!) + +Defines an SCC-split nonlinear system to be solved iteratively. + +## Mathematical Specification of an SCC-Split Nonlinear Problem + +An SCC-Split Nonlinear Problem is a system of nonlinear equations + +```math +f(u,p) = 0 +``` + +with the special property that its Jacobian is in block-lower-triangular +form. In this form, the nonlinear problem can be decomposed into a system +of nonlinear systems. + +```math +f_1(u_1,p) = 0 +f_2(u_2,u_1,p) = 0 +f_3(u_3,u_2,u_1,p) = 0 +\vdots +f_n(u_n,\ldots,u_3,u_2,u_1,p) = 0 +``` + +Splitting the system in this form can have multiple advantages, including: + +* Improved numerical stability and robustness of the solving process +* Improved performance due to using smaller Jacobians + +The SCC-Split Nonlinear Problem is the ordered collection of nonlinear systems +to solve in order solve the system in the optimized split form. + +## Representation + +The representation of the SCCNonlinearProblem is via an ordered collection +of `NonlinearProblem`s, `probs`, with an attached explicit function for pre-processing +a cache. This can be interpreted as follows: + +```math +p_1 = g_1(u,p) +f_1(u_1,p_1) = 0 +p_2 = g_2(u,p) +f_2(u_2,u_1,p_2) = 0 +p_3 = g_3(u,p) +f_3(u_3,u_2,u_1,p_3) = 0 +\vdots +p_n = g_n(u,p) +f_n(u_n,\ldots,u_3,u_2,u_1,p_n) = 0 +``` + +where ``g_i`` is `explicitfuns![i]`. In a computational sense, `explictfuns!` +is instead a mutating function `explictfuns![i](prob.probs[i],sols)` which +updates the values of prob.probs[i] using the previous solutions `sols[i-1]` +and below. + +!!! warn + For the purposes of differentiation, it's assumed that `explictfuns!` does + not modify tunable parameters! + +!!! warn + While `explictfuns![i]` could in theory use `sols[i+1]` in its computation, + these values will not be updated. It is thus the contract of the interface + to not use those values except for as caches to be overriden. + +!!! note + prob.probs[i].p can be aliased with each other as a performance / memory + optimization. If they are aliased before the construction of the `probs` + the runtime guarantees the aliasing behavior is kept. + +## Example + +For the following nonlinear problem: + +```julia +function f(du,u,p) + du[1] = cos(u[2]) - u[1] + du[2] = sin(u[1] + u[2]) + u[2] + du[3] = 2u[4] + u[3] + 1.0 + du[4] = u[5]^2 + u[4] + du[5] = u[3]^2 + u[5] + du[6] = u[1] + u[2] + u[3] + u[4] + u[5] + 2.0u[6] + 2.5u[7] + 1.5u[8] + du[7] = u[1] + u[2] + u[3] + 2.0u[4] + u[5] + 4.0u[6] - 1.5u[7] + 1.5u[8] + du[8] = u[1] + 2.0u[2] + 3.0u[3] + 5.0u[4] + 6.0u[5] + u[6] - u[7] - u[8] +end +prob = NonlinearProblem(f, zeros(8)) +sol = solve(prob) +``` + +The split SCC form is: + +```julia +cache = zeros(3) + +function f1(du,u,cache) + du[1] = cos(u[2]) - u[1] + du[2] = sin(u[1] + u[2]) + u[2] +end +explicitfun1(cache,sols) = nothing +prob1 = NonlinearProblem(NonlinearFunction{true, SciMLBase.NoSpecialize}(f1), zeros(2), cache) +sol1 = solve(prob1, NewtonRaphson()) + +function f2(du,u,cache) + du[1] = 2u[2] + u[1] + 1.0 + du[2] = u[3]^2 + u[2] + du[3] = u[1]^2 + u[3] +end +explicitfun2(cache,sols) = nothing +prob2 = NonlinearProblem(NonlinearFunction{true, SciMLBase.NoSpecialize}(f2), zeros(3), cache) +sol2 = solve(prob2, NewtonRaphson()) + +function f3(du,u,cache) + du[1] = cache[1] + 2.0u[1] + 2.5u[2] + 1.5u[3] + du[2] = cache[2] + 4.0u[1] - 1.5u[2] + 1.5u[3] + du[3] = cache[3] + + u[1] - u[2] - u[3] +end +prob3 = NonlinearProblem(NonlinearFunction{true, SciMLBase.NoSpecialize}(f3), zeros(3), cache) +function explicitfun3(cache,sols) + cache[1] = sols[1][1] + sols[1][2] + sols[2][1] + sols[2][2] + sols[2][3] + cache[2] = sols[1][1] + sols[1][2] + sols[2][1] + 2.0sols[2][2] + sols[2][3] + cache[3] = sols[1][1] + 2.0sols[1][2] + 3.0sols[2][1] + 5.0sols[2][2] + 6.0sols[2][3] +end +explicitfun3(cache,[sol1,sol2]) +sol3 = solve(prob3, NewtonRaphson()) +manualscc = [sol1; sol2; sol3] + +sccprob = SciMLBase.SCCNonlinearProblem([prob1,prob2,prob3], SciMLBase.Void{Any}.([explicitfun1,explicitfun2,explicitfun3])) +``` + +Note that this example aliases the parameters together for a memory-reduced representation. + +## Problem Type + +### Constructors + +### Fields + +* `probs`: the collection of problems to solve +* `explictfuns!`: the explicit functions for mutating the parameter set +""" +mutable struct SCCNonlinearProblem{P, E} + probs::P + explictfuns!::E +end diff --git a/src/solutions/nonlinear_solutions.jl b/src/solutions/nonlinear_solutions.jl index 99acbd8b8..683f625fc 100644 --- a/src/solutions/nonlinear_solutions.jl +++ b/src/solutions/nonlinear_solutions.jl @@ -70,7 +70,7 @@ const SteadyStateSolution = NonlinearSolution get_p(p::AbstractNonlinearSolution) = p.prob.p -function build_solution(prob::AbstractNonlinearProblem, +function build_solution(prob::Union{AbstractNonlinearProblem, SCCNonlinearProblem}, alg, u, resid; calculate_error = true, retcode = ReturnCode.Default, original = nothing,