diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 000000000..8ef51d9d5 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,40 @@ +name: CI +on: + pull_request: + branches: + - master + push: + branches: + - master +jobs: + test: + runs-on: ubuntu-latest + strategy: + matrix: + group: + - Core + - Downstream + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: 1 + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + env: + GROUP: ${{ matrix.group }} + JULIA_NUM_THREADS: 11 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v1 + with: + file: lcov.info diff --git a/.github/workflows/Test.yml b/.github/workflows/Documentation.yml similarity index 53% rename from .github/workflows/Test.yml rename to .github/workflows/Documentation.yml index 1ba42571c..e47f93896 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Documentation.yml @@ -1,4 +1,4 @@ -name: Test +name: Documentation on: push: @@ -9,15 +9,16 @@ on: jobs: build: - runs-on: self-hosted + runs-on: ubuntu-latest steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@latest with: - version: '1.5' + version: '1' - name: Install dependencies - run: julia --project=./ -e 'using Pkg; Pkg.instantiate()' + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - run: julia --project=./ -e 'using Pkg; pkg"test";' + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + run: julia --project=docs/ docs/make.jl diff --git a/README.md b/README.md index be93d8355..d4ddac000 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,20 @@ # NonlinearSolve.jl +[![Github Action CI](https://github.com/SciML/NonlinearSolve.jl/workflows/CI/badge.svg)](https://github.com/SciML/NonlinearSolve.jl/actions) +[![Coverage Status](https://coveralls.io/repos/github/SciML/NonlinearSolve.jl/badge.svg?branch=master)](https://coveralls.io/github/SciML/NonlinearSolve.jl?branch=master) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](http://nlsolve.sciml.ai/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](http://nlsolve.sciml.ai/dev/) +[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac) + Fast implementations of root finding algorithms in Julia that satisfy the SciML common interface. +For information on using the package, +[see the stable documentation](https://mtk.sciml.ai/stable/). Use the +[in-development documentation](https://mtk.sciml.ai/dev/) for the version of +the documentation which contains the unreleased features. + +## High Level Examples + ```julia using NonlinearSolve, StaticArrays @@ -17,37 +30,3 @@ u0 = (1.0, 2.0) # brackets probB = NonlinearProblem(f, u0) sol = solve(probB, Falsi()) ``` - -## Current Algorithms - -### Non-Bracketing - -- `NewtonRaphson()` - -### Bracketing - -- `Falsi()` -- `Bisection()` - -## Features - -Performance is key: the current methods are made to be highly performant on scalar and statically sized small -problems. If you run into any performance issues, please file an issue. - -There is an iterator form of the nonlinear solver which mirrors the DiffEq integrator interface: - -```julia -f(u, p) = u .* u .- 2.0 -u0 = (1.0, 2.0) # brackets -probB = NonlinearProblem(f, u0) -solver = init(probB, Falsi()) # Can iterate the solver object -solver = solve!(solver) -``` - -Note that the `solver` object is actually immutable since we want to make it live on the stack for the sake of performance. - -## Roadmap - -The current algorithms should support automatic differentiation, though improved adjoint overloads are planned -to be added in the current update (which will make use of the `f(u,p)` form). Future updates will include -standard methods for larger scale nonlinear solving like Newton-Krylov methods. diff --git a/docs/.gitignore b/docs/.gitignore new file mode 100644 index 000000000..a303fff20 --- /dev/null +++ b/docs/.gitignore @@ -0,0 +1,2 @@ +build/ +site/ diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 000000000..d93c0ad9a --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,6 @@ +[deps] +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "0.27" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 000000000..b0796ec78 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,31 @@ +using Documenter, NonlinearSolve + +makedocs( + sitename="NonlinearSolve.jl", + authors="Chris Rackauckas", + modules=[NonlinearSolve], + clean=true,doctest=false, + format = Documenter.HTML(#analytics = "UA-90474609-3", + assets = ["assets/favicon.ico"], + canonical="https://nonlinearsolve.sciml.ai/stable/"), + pages=[ + "Home" => "index.md", + "Tutorials" => Any[ + "tutorials/nonlinear.md" + ], + "Basics" => Any[ + "basics/NonlinearProblem.md", + "basics/NonlinearFunctions.md", + "basics/FAQ.md" + ], + "Solvers" => Any[ + "solvers/NonlinearSystemSolvers.md", + "solvers/BracketingSolvers.md" + ] + ] +) + +deploydocs( + repo = "github.com/SciML/NonlinearSolve.jl.git"; + push_preview = true +) diff --git a/docs/src/assets/favicon.ico b/docs/src/assets/favicon.ico new file mode 100644 index 000000000..3c6bd4703 Binary files /dev/null and b/docs/src/assets/favicon.ico differ diff --git a/docs/src/assets/logo.png b/docs/src/assets/logo.png new file mode 100644 index 000000000..6f4c3e261 Binary files /dev/null and b/docs/src/assets/logo.png differ diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md new file mode 100644 index 000000000..c996e4b90 --- /dev/null +++ b/docs/src/basics/FAQ.md @@ -0,0 +1,3 @@ +# Frequently Asked Questions + +Ask more questions. diff --git a/docs/src/basics/NonlinearFunctions.md b/docs/src/basics/NonlinearFunctions.md new file mode 100644 index 000000000..8b2150015 --- /dev/null +++ b/docs/src/basics/NonlinearFunctions.md @@ -0,0 +1,40 @@ +# [NonlinearFunctions and Jacobian Types](@id nonlinearfunctions) + +The SciML ecosystem provides an extensive interface for declaring extra functions +associated with the differential equation's data. In traditional libraries there +is usually only one option: the Jacobian. However, we allow for a large array +of pre-computed functions to speed up the calculations. This is offered via the +`NonlinearFunction` types which can be passed to the problems. + +## Function Type Definitions + +### Function Choice Definitions + +The full interface available to the solvers is as follows: + +- `jac`: The Jacobian of the differential equation with respect to the state + variable `u` at a time `t` with parameters `p`. +- `tgrad`: The gradient of the differential equation with respect to `t` at state + `u` with parameters `p`. +- `paramjac`: The Jacobian of the differential equation with respect to `p` at + state `u` at time `t`. +- `analytic`: Defines an analytical solution using `u0` at time `t` with `p` + which will cause the solvers to return errors. Used for testing. +- `syms`: Allows you to name your variables for automatic names in plots and + other output. + +### NonlinearFunction + +```julia +function NonlinearFunction{iip,true}(f; + analytic=nothing, # (u0,p) + jac=nothing, # (J,u,p) or (u,p) + jvp=nothing, + vjp=nothing, + jac_prototype=nothing, # Type for the Jacobian + sparsity=jac_prototype, + paramjac = nothing, + syms = nothing, + observed = DEFAULT_OBSERVED_NO_TIME, + colorvec = nothing) +``` diff --git a/docs/src/basics/NonlinearProblem.md b/docs/src/basics/NonlinearProblem.md new file mode 100644 index 000000000..5b6b8e72c --- /dev/null +++ b/docs/src/basics/NonlinearProblem.md @@ -0,0 +1,44 @@ +# Nonlinear Problems + +## Mathematical Specification of a Nonlinear Problem + +To define a Nonlinear Problem, you simply need to give the function ``f`` +which defines the nonlinear system: + +```math +f(u,p) = 0 +``` + +and an initial guess ``u₀`` of where `f(u,p)=0`. `f` should be specified as `f(u,p)` +(or in-place as `f(du,u,p)`), and `u₀` should be an AbstractArray (or number) +whose geometry matches the desired geometry of `u`. Note that we are not limited +to numbers or vectors for `u₀`; one is allowed to provide `u₀` as arbitrary +matrices / higher dimension tensors as well. + +## Problem Type + +### Constructors + +```julia +NonlinearProblem(f::NonlinearFunction,u0,p=NullParameters();kwargs...) +NonlinearProblem{isinplace}(f,u0,p=NullParameters();kwargs...) +``` + +`isinplace` optionally sets whether the function is inplace or not. This is +determined automatically, but not inferred. + +Parameters are optional, and if not given then a `NullParameters()` singleton +will be used which will throw nice errors if you try to index non-existent +parameters. Any extra keyword arguments are passed on to the solvers. For example, +if you set a `callback` in the problem, then that `callback` will be added in +every solve call. + +For specifying Jacobians and mass matrices, see the [NonlinearFunctions](@ref nonlinearfunctions) +page. + +### Fields + +* `f`: The function in the ODE. +* `u0`: The initial guess for the steady state. +* `p`: The parameters for the problem. Defaults to `NullParameters` +* `kwargs`: The keyword arguments passed onto the solves. diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 000000000..b38e1ae53 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,42 @@ +# NonlinearSolve.jl: High-Performance Unified Nonlinear Solvers + +NonlinaerSolve.jl is a unified interface for the nonlinear solving packages of +Julia. It includes its own high-performance nonlinear solves which include the +ability to swap out to fast direct and iterative linear solvers, along with the +ability to use sparse automatic differentiation for Jacobian construction and +Jacobian-vector products. It interfaces with other packages of the Julia ecosystem +to make it easy to test alternative solver packages and pass small types to +control algorithm swapping. It interfaces with the +[ModelingToolkit.jl](https://mtk.sciml.ai/dev/) world of symbolic modeling to +allow for automatically generating high performance code. + +Performance is key: the current methods are made to be highly performant on +scalar and statically sized small problems, with options for large-scale systems. +If you run into any performance issues, please file an issue. + +## Installation + +To install NonlinearSolve.jl, use the Julia package manager: + +```julia +using Pkg +Pkg.add("NonlinearSolve") +``` + +## Contributing + +- Please refer to the + [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md) + for guidance on PRs, issues, and other matters relating to contributing to ModelingToolkit. +- There are a few community forums: + - The #diffeq-bridged channel in the [Julia Slack](https://julialang.org/slack/) + - [JuliaDiffEq](https://gitter.im/JuliaDiffEq/Lobby) on Gitter + - On the Julia Discourse forums (look for the [modelingtoolkit tag](https://discourse.julialang.org/tag/modelingtoolkit) + - See also [SciML Community page](https://sciml.ai/community/) + +## Roadmap + +The current algorithms should support automatic differentiation, though improved +adjoint overloads are planned to be added in the current update (which will make +use of the `f(u,p)` form). Future updates will include standard methods for +larger scale nonlinear solving like Newton-Krylov methods. diff --git a/docs/src/solvers/BracketingSolvers.md b/docs/src/solvers/BracketingSolvers.md new file mode 100644 index 000000000..293ddc28b --- /dev/null +++ b/docs/src/solvers/BracketingSolvers.md @@ -0,0 +1,20 @@ +# Bracketing Solvers + +`solve(prob::NonlinearProblem,alg;kwargs)` + +Solves for ``f(u)=0`` in the problem defined by `prob` using the algorithm +`alg`. If no algorithm is given, a default algorithm will be chosen. + +This page is solely focused on the bracketing methods for scalar nonlinear equations. + +## Recommended Methods + +`Falsi()` can have a faster convergence and is discretely differentiable, but is +less stable than `Bisection`. + +## Full List of Methods + +### NonlinearSolve.jl + +- `Falsi` : A non-allocating regula falsi method +- `Bisection`: A common bisection method diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md new file mode 100644 index 000000000..17669625b --- /dev/null +++ b/docs/src/solvers/NonlinearSystemSolvers.md @@ -0,0 +1,117 @@ +# Nonlinear System Solvers + +`solve(prob::NonlinearProblem,alg;kwargs)` + +Solves for ``f(u)=0`` in the problem defined by `prob` using the algorithm +`alg`. If no algorithm is given, a default algorithm will be chosen. + +This page is solely focused on the methods for nonlinear systems. + +## Recommended Methods + +`NewtonRaphson` is a good choice for most problems. It is non-allocating on +static arrays and thus really well-optimized for small systems, while for large +systems it can make use of sparsity patterns for sparse automatic differentiation +and sparse linear solving of very large systems. That said, as a classic Newton +method, it's stability region can be smaller than other methods. `NLSolveJL`'s +`:trust_region` method can be a good choice for high stability, along with +`CMINPACK`. + +For a system which is very non-stiff (i.e. the condition number of the Jacobian +is small, or the eigenvalues of the Jacobian are within a few orders of magnitude), +then `NLSolveJL`'s `:anderson` can be a good choice. + +## Full List of Methods + +### NonlinearSolve.jl + +These are the core solvers. + +- `NewtonRaphson(;autodiff=true,chunk_size=12,diff_type=Val{:forward},linsolve=DEFAULT_LINSOLVE)`: + A Newton-Raphson method with swappable nonlinear solvers and autodiff methods + for high performance on large and sparse systems. When used on objects like + static arrays, this method is non-allocating. + +### SciMLNLSolve.jl + +This is a wrapper package to import solvers from other packages into this interface. +Note that these solvers do not come by default, and thus one needs to install +the package before using these solvers: + +```julia +]add SciMLNLSolve +using SciMLNLSolve +``` + +- `CMINPACK()`: A wrapper for using the classic MINPACK method through [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) +- `NLSolveJL()`: A wrapper for [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) + +```julia +NLSolveJL(; + method=:trust_region, + autodiff=:central, + store_trace=false, + extended_trace=false, + linesearch=LineSearches.Static(), + linsolve=(x, A, b) -> copyto!(x, A\b), + factor = one(Float64), + autoscale=true, + m=10, + beta=one(Float64), + show_trace=false, + ) +``` + +Choices for methods in `NLSolveJL`: + +- `:fixedpoint`: fixed point iteration +- `:anderson`: Anderson accelerated fixed point iteration +- `:newton`: Classical Newton method with an optional line search +- `:trust_region`: Trust region Newton method. The default + +For more information on these arguments, consult the +[NLsolve.jl documentation](https://github.com/JuliaNLSolvers/NLsolve.jl) + +### Sundials.jl + +This is a wrapper package for the SUNDIALS C library, specifically the KINSOL +nonlinear solver included in that ecosystem. Note that these solvers do not come +by default, and thus one needs to install the package before using these solvers: + +```julia +]add Sundials +using Sundials +``` + +- `KINSOL`: The KINSOL method of the SUNDIALS C library + +```julia +KINSOL(; + linear_solver = :Dense, + jac_upper = 0, + jac_lower = 0, + userdata = nothing, +) +``` + +The choices for the linear solver are: + +- `:Dense` - A dense linear solver. +- `:Band` - A solver specialized for banded Jacobians. If used, you must set the + position of the upper and lower non-zero diagonals via `jac_upper` and + `jac_lower`. +- `:LapackDense` - A version of the dense linear solver that uses the Julia-provided + OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than + `:Dense` on larger systems but has noticable overhead on smaller (<100 ODE) systems. +- `:LapackBand` - A version of the banded linear solver that uses the Julia-provided + OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than + `:Band` on larger systems but has noticable overhead on smaller (<100 ODE) systems. +- `:Diagonal` - This method is specialized for diagonal Jacobians. +- `:GMRES` - A GMRES method. Recommended first choice Krylov method +- `:BCG` - A Biconjugate gradient method. +- `:PCG` - A preconditioned conjugate gradient method. Only for symmetric + linear systems. +- `:TFQMR` - A TFQMR method. +- `:KLU` - A sparse factorization method. Requires that the user specifies a + Jacobian. The Jacobian must be set as a sparse matrix in the `ODEProblem` + type. diff --git a/docs/src/tutorials/iterator_interface.md b/docs/src/tutorials/iterator_interface.md new file mode 100644 index 000000000..f66c6c768 --- /dev/null +++ b/docs/src/tutorials/iterator_interface.md @@ -0,0 +1,14 @@ +# Nonlinear Solver Iterator Interface + +There is an iterator form of the nonlinear solver which mirrors the DiffEq integrator interface: + +```julia +f(u, p) = u .* u .- 2.0 +u0 = (1.0, 2.0) # brackets +probB = NonlinearProblem(f, u0) +solver = init(probB, Falsi()) # Can iterate the solver object +solver = solve!(solver) +``` + +Note that the `solver` object is actually immutable since we want to make it +live on the stack for the sake of performance. diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md new file mode 100644 index 000000000..3974ecfc1 --- /dev/null +++ b/docs/src/tutorials/nonlinear.md @@ -0,0 +1,32 @@ +# Solving Nonlinear Systems + +A nonlinear system $$f(u) = 0$$ is specified by defining a function `f(u,p)`, +where `p` are the parameters of the system. For example, the following solves +the vector equation $$f(u) = u^2 - p$$ for a vector of equations: + +```julia +using NonlinearSolve, StaticArrays + +f(u,p) = u .* u .- p +u0 = @SVector[1.0, 1.0] +p = 2.0 +probN = NonlinearProblem{false}(f, u0, p) +solver = solve(probN, NewtonRaphson(), tol = 1e-9) +``` + +where `u0` is the initial condition for the rootfind. Native NonlinearSolve.jl +solvers use the given type of `u0` to determine the type used within the solver +and the return. Note that the parameters `p` can be any type, but most be an +AbstractArray for automatic differentiation. + +## Using Bracketing Methods + +For scalar rootfinding problems, bracketing methods exist. In this case one passes +a bracket instead of an initial condition, for example: + +```julia +f(u, p) = u .* u .- 2.0 +u0 = (1.0, 2.0) # brackets +probB = NonlinearProblem(f, u0) +sol = solve(probB, Falsi()) +```