From 6c0e4e1dd139f51bc6f6e4f0a93c5560d6cafa16 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:44:15 +0200 Subject: [PATCH 01/30] Implement status and params function for nonlinear solvers. --- src/solvers/nonlinear/abstract_newton_solver.jl | 2 ++ src/solvers/nonlinear/nonlinear_solvers.jl | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/solvers/nonlinear/abstract_newton_solver.jl b/src/solvers/nonlinear/abstract_newton_solver.jl index 09c0972d9..3901329fc 100644 --- a/src/solvers/nonlinear/abstract_newton_solver.jl +++ b/src/solvers/nonlinear/abstract_newton_solver.jl @@ -29,6 +29,8 @@ end function computeJacobian(s::AbstractNewtonSolver) computeJacobian(s.x, s.J, s.Jparams) end +status(solver::AbstractNewtonSolver) = solver.status +params(solver::AbstractNewtonSolver) = solver.params function check_jacobian(s::AbstractNewtonSolver) println("Condition Number of Jacobian: ", cond(s.J)) diff --git a/src/solvers/nonlinear/nonlinear_solvers.jl b/src/solvers/nonlinear/nonlinear_solvers.jl index 078115cde..9cd7d233e 100644 --- a/src/solvers/nonlinear/nonlinear_solvers.jl +++ b/src/solvers/nonlinear/nonlinear_solvers.jl @@ -4,6 +4,8 @@ using Printf abstract type NonlinearSolver{T} end solve!(s::NonlinearSolver) = error("solve! not implemented for $(typeof(s))") +status(s::NonlinearSolver) = error("status not implemented for $(typeof(s))") +params(s::NonlinearSolver) = error("params not implemented for $(typeof(s))") function solve!(s::NonlinearSolver{T}, x₀::Vector{T}) where {T} setInitialConditions!(s, x₀) From 429778a4bc9ad6a6ec6aa3dc959b5b065ae6572a Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:48:49 +0200 Subject: [PATCH 02/30] Implement get_solver_status functions for nonlinear solvers returning a dictionary with number of iterations, residuals, etc. (+1 squashed commit) Squashed commits: [8f605d5] export solver status functions. --- src/Solvers.jl | 1 + src/solvers/nonlinear/nonlinear_solvers.jl | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/src/Solvers.jl b/src/Solvers.jl index bd9afce93..ef877bca9 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -23,6 +23,7 @@ module Solvers NewtonSolver, QuasiNewtonSolver, residual_initial!, residual_absolute!, residual_relative!, print_solver_status, check_solver_converged, check_solver_status, + get_solver_status, get_solver_status!, solve! include("solvers/nonlinear/nonlinear_solvers.jl") diff --git a/src/solvers/nonlinear/nonlinear_solvers.jl b/src/solvers/nonlinear/nonlinear_solvers.jl index 9cd7d233e..d22007751 100644 --- a/src/solvers/nonlinear/nonlinear_solvers.jl +++ b/src/solvers/nonlinear/nonlinear_solvers.jl @@ -113,6 +113,26 @@ function check_solver_status(status::NonlinearSolverStatus, params::NonlinearSol end end +function get_solver_status!(status::NonlinearSolverStatus, params::NonlinearSolverParameters, status_dict::Dict) + status_dict[:nls_niter] = status.i + status_dict[:nls_atol] = status.rₐ + status_dict[:nls_rtol] = status.rᵣ + status_dict[:nls_stol] = status.rₛ + status_dict[:nls_converged] = check_solver_converged(status, params) + return status_dict +end + +get_solver_status!(solver::NonlinearSolver{T}, status_dict::Dict) where {T} = + get_solver_status!(status(solver), params(solver), status_dict) + +get_solver_status(solver::NonlinearSolver{T}) where {T} = get_solver_status!(solver, + Dict(:nls_niter => 0, + :nls_atol => zero(T), + :nls_rtol => zero(T), + :nls_stol => zero(T), + :nls_converged => false) + ) + function getLinearSolver(x::AbstractVector{T}) where {T} n = length(x) From f0d2823fc301a515c432f96d16beee22f9fdecb9 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:45:17 +0200 Subject: [PATCH 03/30] Minor cleanup in abstract Newton solver. --- src/solvers/nonlinear/abstract_newton_solver.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/solvers/nonlinear/abstract_newton_solver.jl b/src/solvers/nonlinear/abstract_newton_solver.jl index 3901329fc..b850a3397 100644 --- a/src/solvers/nonlinear/abstract_newton_solver.jl +++ b/src/solvers/nonlinear/abstract_newton_solver.jl @@ -21,17 +21,15 @@ abstract type AbstractNewtonSolver{T} <: NonlinearSolver{T} end status::NonlinearSolverStatus{T} end - function setInitialConditions!(s::AbstractNewtonSolver{T}, x₀::Vector{T}) where {T} s.x[:] = x₀ end -function computeJacobian(s::AbstractNewtonSolver) - computeJacobian(s.x, s.J, s.Jparams) -end status(solver::AbstractNewtonSolver) = solver.status params(solver::AbstractNewtonSolver) = solver.params +computeJacobian(s::AbstractNewtonSolver) = computeJacobian(s.x, s.J, s.Jparams) + function check_jacobian(s::AbstractNewtonSolver) println("Condition Number of Jacobian: ", cond(s.J)) println("Determinant of Jacobian: ", det(s.J)) From f520bceb789c2c5b5a589c0eac532f5fa46be2e9 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:46:10 +0200 Subject: [PATCH 04/30] Cleanup atomic solution constructor calls. --- src/solutions/atomic_solution.jl | 12 ++++++------ src/solutions/solution_dae.jl | 2 +- src/solutions/solution_ode.jl | 2 +- src/solutions/solution_pdae.jl | 2 +- src/solutions/solution_pode.jl | 2 +- src/solutions/solution_psde.jl | 2 +- src/solutions/solution_sde.jl | 2 +- 7 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/solutions/atomic_solution.jl b/src/solutions/atomic_solution.jl index dc1d899a8..3ca36db3c 100644 --- a/src/solutions/atomic_solution.jl +++ b/src/solutions/atomic_solution.jl @@ -7,32 +7,32 @@ abstract type AtomicSolution{dType <: Number, tType <: Real} end "Create AtomicSolution for ODE." function AtomicSolution(equation::AbstractEquationODE{DT,TT}) where {DT,TT} - AtomicSolutionODE{DT,TT}(ndims(equation)) + AtomicSolutionODE(DT, TT, ndims(equation)) end "Create AtomicSolution for partitioned ODE." function AtomicSolution(equation::AbstractEquationPODE{DT,TT}) where {DT,TT} - AtomicSolutionPODE{DT,TT}(ndims(equation)) + AtomicSolutionPODE(DT, TT, ndims(equation)) end "Create AtomicSolution for DAE." function AtomicSolution(equation::AbstractEquationDAE{DT,TT}) where {DT,TT} - AtomicSolutionDAE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionDAE(DT, TT, ndims(equation), equation.m) end "Create AtomicSolution for partitioned DAE." function AtomicSolution(equation::AbstractEquationPDAE{DT,TT}) where {DT,TT} - AtomicSolutionPDAE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionPDAE(DT, TT, ndims(equation), equation.m) end "Create AtomicSolution for SDE." function AtomicSolution(equation::AbstractEquationSDE{DT,TT}) where {DT,TT} - AtomicSolutionSDE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionSDE(DT, TT, ndims(equation), equation.m) end "Create AtomicSolution for PSDE." function AtomicSolution(equation::AbstractEquationPSDE{DT,TT}) where {DT,TT} - AtomicSolutionPSDE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionPSDE(DT, TT, ndims(equation), equation.m) end "Print error for AtomicSolutions of equations not implemented, yet." diff --git a/src/solutions/solution_dae.jl b/src/solutions/solution_dae.jl index ed9bb8d1a..b1515cc40 100644 --- a/src/solutions/solution_dae.jl +++ b/src/solutions/solution_dae.jl @@ -157,7 +157,7 @@ Base.:(==)(sol1::SolutionDAE{DT1,TT1,N1}, sol2::SolutionDAE{DT2,TT2,N2}) where { "Create AtomicSolution for DAE." function AtomicSolution(solution::SolutionDAE{DT,TT}) where {DT,TT} - AtomicSolutionDAE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionDAE(DT, TT, solution.nd, solution.nm) end diff --git a/src/solutions/solution_ode.jl b/src/solutions/solution_ode.jl index 513851aa1..b94556d17 100644 --- a/src/solutions/solution_ode.jl +++ b/src/solutions/solution_ode.jl @@ -142,7 +142,7 @@ Base.:(==)(sol1::SolutionODE{DT1,TT1,N1}, sol2::SolutionODE{DT2,TT2,N2}) where { "Create AtomicSolution for ODE." function AtomicSolution(solution::SolutionODE{DT,TT}) where {DT,TT} - AtomicSolutionODE{DT,TT}(solution.nd) + AtomicSolutionODE(DT, TT, solution.nd) end diff --git a/src/solutions/solution_pdae.jl b/src/solutions/solution_pdae.jl index 09a669837..e99c254cc 100644 --- a/src/solutions/solution_pdae.jl +++ b/src/solutions/solution_pdae.jl @@ -163,7 +163,7 @@ Base.:(==)(sol1::SolutionPDAE{DT1,TT1,N1}, sol2::SolutionPDAE{DT2,TT2,N2}) where "Create AtomicSolution for partitioned DAE." function AtomicSolution(solution::SolutionPDAE{DT,TT}) where {DT,TT} - AtomicSolutionPDAE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionPDAE(DT, TT, solution.nd, solution.nm) end diff --git a/src/solutions/solution_pode.jl b/src/solutions/solution_pode.jl index 0902b200a..375100239 100644 --- a/src/solutions/solution_pode.jl +++ b/src/solutions/solution_pode.jl @@ -147,7 +147,7 @@ Base.:(==)(sol1::SolutionPODE, sol2::SolutionPODE) = ( "Create AtomicSolution for partitioned ODE." function AtomicSolution(solution::SolutionPODE{DT,TT}) where {DT,TT} - AtomicSolutionPODE{DT,TT}(solution.nd) + AtomicSolutionPODE(DT, TT, solution.nd) end diff --git a/src/solutions/solution_psde.jl b/src/solutions/solution_psde.jl index e438c4b60..2b37a5946 100644 --- a/src/solutions/solution_psde.jl +++ b/src/solutions/solution_psde.jl @@ -259,7 +259,7 @@ Base.:(==)(sol1::SolutionPSDE{DT1,TT1,NQ1,NW1,C1}, sol2::SolutionPSDE{DT2,TT2,NQ "Create AtomicSolution for PSDE." function AtomicSolution(solution::SolutionPSDE{DT,TT}) where {DT,TT} - AtomicSolutionPSDE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionPSDE(DT, TT, solution.nd, solution.nm) end diff --git a/src/solutions/solution_sde.jl b/src/solutions/solution_sde.jl index 91060d0c1..5f19ba4d0 100644 --- a/src/solutions/solution_sde.jl +++ b/src/solutions/solution_sde.jl @@ -250,7 +250,7 @@ Base.:(==)(sol1::SolutionSDE{DT1,TT1,NQ1,NW1,C1}, sol2::SolutionSDE{DT2,TT2,NQ2, "Create AtomicSolution for SDE." function AtomicSolution(solution::SolutionSDE{DT,TT}) where {DT,TT} - AtomicSolutionSDE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionSDE(DT, TT, solution.nd, solution.nm) end From 51c8164e24ea9f9e7ff41fa269d53d37be2d1066 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:49:39 +0200 Subject: [PATCH 05/30] Implement nconstraints functions for differential algebraic equations and corresponding integrators. Clean up ndims function. (+2 squashed commits) Squashed commits: [e9863cb] Add ndims functions for Galerkin and splitting integrators. [9d4140b] Add fallback nconstraints function for integrators. --- src/CommonFunctions.jl | 3 +++ src/equations/hdae.jl | 4 ++-- src/equations/idae.jl | 4 ++-- src/equations/pdae.jl | 4 ++-- src/equations/spdae.jl | 4 ++-- src/equations/vdae.jl | 4 ++-- src/integrators/abstract_integrator.jl | 1 + src/integrators/cgvi/integrators_cgvi.jl | 2 ++ src/integrators/dgvi/integrators_dgvi.jl | 2 ++ src/integrators/dgvi/integrators_dgvi_experimental.jl | 1 + src/integrators/dgvi/integrators_dgvi_path_integral.jl | 1 + src/integrators/dgvi/integrators_dgvi_projection_final.jl | 1 + src/integrators/dgvi/integrators_dgvi_projection_initial.jl | 1 + src/integrators/rk/integrators_dirk.jl | 2 +- src/integrators/rk/integrators_eprk.jl | 2 +- src/integrators/rk/integrators_erk.jl | 2 +- src/integrators/rk/integrators_firk.jl | 2 +- src/integrators/rk/integrators_flrk.jl | 2 +- src/integrators/rk/integrators_iprk.jl | 2 +- src/integrators/rk/integrators_pglrk.jl | 4 ++-- src/integrators/spark/integrators_hpark.jl | 3 +++ src/integrators/spark/integrators_hspark.jl | 3 +++ src/integrators/spark/integrators_hspark_primary.jl | 3 +++ src/integrators/spark/integrators_hspark_secondary.jl | 3 +++ src/integrators/spark/integrators_slrk.jl | 3 +++ src/integrators/spark/integrators_vpark.jl | 3 +++ src/integrators/spark/integrators_vspark.jl | 3 +++ src/integrators/spark/integrators_vspark_primary.jl | 3 +++ src/integrators/spark/integrators_vspark_secondary.jl | 3 +++ src/integrators/splitting/integrators_composition.jl | 1 + src/integrators/splitting/integrators_exact_ode.jl | 1 + src/integrators/splitting/integrators_splitting.jl | 1 + 32 files changed, 60 insertions(+), 18 deletions(-) diff --git a/src/CommonFunctions.jl b/src/CommonFunctions.jl index c2d865b01..6e979db7e 100644 --- a/src/CommonFunctions.jl +++ b/src/CommonFunctions.jl @@ -33,4 +33,7 @@ module CommonFunctions eachsample() = nothing eachtimestep() = nothing + export nconstraints + nconstraints() = nothing + end diff --git a/src/equations/hdae.jl b/src/equations/hdae.jl index 1397331ad..271edd755 100644 --- a/src/equations/hdae.jl +++ b/src/equations/hdae.jl @@ -123,8 +123,8 @@ function Base.similar(dae::HDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra parameters=parameters, periodicity=periodicity) end -@inline Base.ndims(dae::HDAE) = dae.d - +@inline Base.ndims(equation::HDAE) = equation.d +@inline CommonFunctions.nconstraints(equation::HDAE) = equation.m @inline CommonFunctions.periodicity(equation::HDAE) = equation.periodicity function get_function_tuple(equation::HDAE{DT,TT,VT,FT,UT,GT,U̅T,G̅T,ϕT,ψT,HT,Nothing}) where {DT, TT, VT, FT, UT, GT, U̅T, G̅T, ϕT, ψT, HT} diff --git a/src/equations/idae.jl b/src/equations/idae.jl index de22bbbca..c1dee93b8 100644 --- a/src/equations/idae.jl +++ b/src/equations/idae.jl @@ -122,8 +122,8 @@ function Base.similar(dae::IDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra IDAE(dae.ϑ, dae.f, dae.u, dae.g, dae.ϕ, t₀, q₀, p₀, λ₀; h=h, v=v, parameters=parameters, periodicity=periodicity) end -@inline Base.ndims(dae::IDAE) = dae.d - +@inline Base.ndims(equation::IDAE) = equation.d +@inline CommonFunctions.nconstraints(equation::IDAE) = equation.m @inline CommonFunctions.periodicity(equation::IDAE) = equation.periodicity function get_function_tuple(equation::IDAE{DT,TT,ϑT,FT,UT,GT,ϕT,HT,VT,Nothing}) where {DT, TT, ϑT, FT, UT, GT, ϕT, HT, VT} diff --git a/src/equations/pdae.jl b/src/equations/pdae.jl index 70b467091..ba60b57fb 100644 --- a/src/equations/pdae.jl +++ b/src/equations/pdae.jl @@ -117,8 +117,8 @@ function Base.similar(dae::PDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra PDAE(dae.v, dae.f, dae.u, dae.g, dae.ϕ, t₀, q₀, p₀, λ₀; h=h, parameters=parameters, periodicity=periodicity) end -@inline Base.ndims(dae::PDAE) = dae.d - +@inline Base.ndims(equation::PDAE) = equation.d +@inline CommonFunctions.nconstraints(equation::PDAE) = equation.m @inline CommonFunctions.periodicity(equation::PDAE) = equation.periodicity function get_function_tuple(equation::PDAE{DT,TT,VT,FT,UT,GT,ϕT,HT,Nothing}) where {DT, TT, VT, FT, UT, GT, ϕT, HT} diff --git a/src/equations/spdae.jl b/src/equations/spdae.jl index 726bde7fb..7efd95aba 100644 --- a/src/equations/spdae.jl +++ b/src/equations/spdae.jl @@ -105,8 +105,8 @@ function Base.similar(dae::SPDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstr parameters=parameters, periodicity=periodicity) end -@inline Base.ndims(dae::SPDAE) = dae.d - +@inline Base.ndims(equation::SPDAE) = equation.d +@inline CommonFunctions.nconstraints(equation::SPDAE) = equation.m @inline CommonFunctions.periodicity(equation::SPDAE) = equation.periodicity function get_function_tuple(equation::SPDAE{DT,TT,VT,FT,ϕT,ψT,Nothing}) where {DT, TT, VT, FT, ϕT, ψT} diff --git a/src/equations/vdae.jl b/src/equations/vdae.jl index 034428208..2ab8fdbe3 100644 --- a/src/equations/vdae.jl +++ b/src/equations/vdae.jl @@ -175,8 +175,8 @@ function Base.similar(dae::VDAE, t₀::TT, q₀::AbstractArray{DT}, p₀::Abstra h=h, v=v, Ω=Ω, ∇H=∇H, parameters=parameters, periodicity=periodicity) end -@inline Base.ndims(dae::VDAE) = dae.d - +@inline Base.ndims(equation::VDAE) = equation.d +@inline CommonFunctions.nconstraints(equation::VDAE) = equation.m @inline CommonFunctions.periodicity(equation::VDAE) = equation.periodicity diff --git a/src/integrators/abstract_integrator.jl b/src/integrators/abstract_integrator.jl index e50eb9cd1..60a52a039 100644 --- a/src/integrators/abstract_integrator.jl +++ b/src/integrators/abstract_integrator.jl @@ -23,6 +23,7 @@ abstract type SPSDEIntegrator{dType, tType} <: StochasticIntegrator{dType, tType equation(integrator::Integrator) = error("equation() not implemented for ", typeof(integrator)) timestep(integrator::Integrator) = error("timestep() not implemented for ", typeof(integrator)) Base.ndims(integrator::Integrator) = error("ndims() not implemented for ", typeof(integrator)) +CommonFunctions.nconstraints(integrator::Integrator) = error("nconstraints() not implemented for ", typeof(integrator)) eachdim(integrator::Integrator) = 1:ndims(integrator) diff --git a/src/integrators/cgvi/integrators_cgvi.jl b/src/integrators/cgvi/integrators_cgvi.jl index c100a8b28..93a6649c3 100644 --- a/src/integrators/cgvi/integrators_cgvi.jl +++ b/src/integrators/cgvi/integrators_cgvi.jl @@ -171,6 +171,8 @@ end @inline equation(integrator::IntegratorCGVI, i::Symbol) = integrator.params.equs[i] @inline equations(integrator::IntegratorCGVI) = integrator.params.equs @inline timestep(integrator::IntegratorCGVI) = integrator.params.Δt +@inline Base.ndims(::IntegratorCGVI{DT,TT,D}) where {DT,TT,D} = D + function IntegratorCache(int::IntegratorCGVI{DT,TT}) where {DT,TT} IntegratorCacheCGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}() diff --git a/src/integrators/dgvi/integrators_dgvi.jl b/src/integrators/dgvi/integrators_dgvi.jl index b5d644e5c..85f82b008 100644 --- a/src/integrators/dgvi/integrators_dgvi.jl +++ b/src/integrators/dgvi/integrators_dgvi.jl @@ -424,6 +424,8 @@ end @inline equation(integrator::IntegratorDGVI, i::Symbol) = integrator.params.equs[i] @inline equations(integrator::IntegratorDGVI) = integrator.params.equs @inline timestep(integrator::IntegratorDGVI) = integrator.params.Δt +@inline Base.ndims(::IntegratorDGVI{DT,TT,D}) where {DT,TT,D} = D + function IntegratorCache(int::IntegratorDGVI{DT,TT}) where {DT,TT} IntegratorCacheDGVI{DT, TT, ndims(int), nbasis(int.basis), nnodes(int.quadrature)}() diff --git a/src/integrators/dgvi/integrators_dgvi_experimental.jl b/src/integrators/dgvi/integrators_dgvi_experimental.jl index 3c4636ce6..cf5813cb9 100644 --- a/src/integrators/dgvi/integrators_dgvi_experimental.jl +++ b/src/integrators/dgvi/integrators_dgvi_experimental.jl @@ -177,6 +177,7 @@ end @inline equation(integrator::IntegratorDGVIEXP, i::Symbol) = integrator.params.equs[i] @inline equations(integrator::IntegratorDGVIEXP) = integrator.params.equs @inline timestep(integrator::IntegratorDGVIEXP) = integrator.params.Δt +@inline Base.ndims(::IntegratorDGVIEXP{DT,TT,D}) where {DT,TT,D} = D function update_params!(params::ParametersDGVIEXP, int::IntegratorDGVIEXP) diff --git a/src/integrators/dgvi/integrators_dgvi_path_integral.jl b/src/integrators/dgvi/integrators_dgvi_path_integral.jl index 0f3f8dde7..7659c3c36 100644 --- a/src/integrators/dgvi/integrators_dgvi_path_integral.jl +++ b/src/integrators/dgvi/integrators_dgvi_path_integral.jl @@ -321,6 +321,7 @@ end @inline equation(integrator::IntegratorDGVIPI, i::Symbol) = integrator.params.equs[i] @inline equations(integrator::IntegratorDGVIPI) = integrator.params.equs @inline timestep(integrator::IntegratorDGVIPI) = integrator.params.Δt +@inline Base.ndims(::IntegratorDGVIPI{DT,TT,D}) where {DT,TT,D} = D function update_params!(params::ParametersDGVIPI, int::IntegratorDGVIPI) diff --git a/src/integrators/dgvi/integrators_dgvi_projection_final.jl b/src/integrators/dgvi/integrators_dgvi_projection_final.jl index 034098810..fdf665301 100644 --- a/src/integrators/dgvi/integrators_dgvi_projection_final.jl +++ b/src/integrators/dgvi/integrators_dgvi_projection_final.jl @@ -177,6 +177,7 @@ end @inline equation(integrator::IntegratorDGVIP1, i::Symbol) = integrator.params.equs[i] @inline equations(integrator::IntegratorDGVIP1) = integrator.params.equs @inline timestep(integrator::IntegratorDGVIP1) = integrator.params.Δt +@inline Base.ndims(::IntegratorDGVIP1{DT,TT,D}) where {DT,TT,D} = D function update_params!(params::ParametersDGVIP1, int::IntegratorDGVIP1) diff --git a/src/integrators/dgvi/integrators_dgvi_projection_initial.jl b/src/integrators/dgvi/integrators_dgvi_projection_initial.jl index a08634efd..7bea8d150 100644 --- a/src/integrators/dgvi/integrators_dgvi_projection_initial.jl +++ b/src/integrators/dgvi/integrators_dgvi_projection_initial.jl @@ -188,6 +188,7 @@ end @inline equation(integrator::IntegratorDGVIP0, i::Symbol) = integrator.params.equs[i] @inline equations(integrator::IntegratorDGVIP0) = integrator.params.equs @inline timestep(integrator::IntegratorDGVIP0) = integrator.params.Δt +@inline Base.ndims(::IntegratorDGVIP0{DT,TT,D}) where {DT,TT,D} = D function update_params!(params::ParametersDGVIP0, int::IntegratorDGVIP0) diff --git a/src/integrators/rk/integrators_dirk.jl b/src/integrators/rk/integrators_dirk.jl index 4dd562b26..71a39d1ed 100644 --- a/src/integrators/rk/integrators_dirk.jl +++ b/src/integrators/rk/integrators_dirk.jl @@ -119,7 +119,7 @@ struct IntegratorDIRK{DT, TT, D, S, PT <: ParametersDIRK{DT,TT}, end -@inline Base.ndims(int::IntegratorDIRK{DT,TT,D,S}) where {DT,TT,D,S} = D +@inline Base.ndims(::IntegratorDIRK{DT,TT,D,S}) where {DT,TT,D,S} = D @inline has_initial_guess(int::IntegratorDIRK) = true diff --git a/src/integrators/rk/integrators_eprk.jl b/src/integrators/rk/integrators_eprk.jl index 239e54409..4d95c707e 100644 --- a/src/integrators/rk/integrators_eprk.jl +++ b/src/integrators/rk/integrators_eprk.jl @@ -122,7 +122,7 @@ struct IntegratorEPRK{DT, TT, D, S, ET <: NamedTuple} <: IntegratorPRK{DT,TT} end -@inline Base.ndims(int::IntegratorEPRK{DT,TT,D,S}) where {DT,TT,D,S} = D +@inline Base.ndims(::IntegratorEPRK{DT,TT,D,S}) where {DT,TT,D,S} = D "Compute Q stages of explicit partitioned Runge-Kutta methods." diff --git a/src/integrators/rk/integrators_erk.jl b/src/integrators/rk/integrators_erk.jl index d415e0f09..8773b955d 100644 --- a/src/integrators/rk/integrators_erk.jl +++ b/src/integrators/rk/integrators_erk.jl @@ -117,7 +117,7 @@ struct IntegratorERK{DT, TT, D, S, ET} <: IntegratorRK{DT,TT} end -@inline Base.ndims(int::IntegratorERK{DT,TT,PT,D,S}) where {DT,TT,PT,D,S} = D +@inline Base.ndims(::IntegratorERK{DT,TT,D,S}) where {DT,TT,D,S} = D "Integrate ODE with explicit Runge-Kutta integrator." diff --git a/src/integrators/rk/integrators_firk.jl b/src/integrators/rk/integrators_firk.jl index 211554d62..1d806b8f6 100644 --- a/src/integrators/rk/integrators_firk.jl +++ b/src/integrators/rk/integrators_firk.jl @@ -143,7 +143,7 @@ struct IntegratorFIRK{DT, TT, D, S, PT <: ParametersFIRK{DT,TT}, end -@inline Base.ndims(int::IntegratorFIRK{DT,TT,D,S}) where {DT,TT,D,S} = D +@inline Base.ndims(::IntegratorFIRK{DT,TT,D,S}) where {DT,TT,D,S} = D function initialize!(int::IntegratorFIRK, sol::AtomicSolutionODE) diff --git a/src/integrators/rk/integrators_flrk.jl b/src/integrators/rk/integrators_flrk.jl index e8fa9fe34..47aa272a5 100644 --- a/src/integrators/rk/integrators_flrk.jl +++ b/src/integrators/rk/integrators_flrk.jl @@ -103,7 +103,7 @@ struct IntegratorFLRK{DT, TT, D, S, PT <: ParametersFLRK{DT,TT}, end -@inline Base.ndims(int::IntegratorFLRK{DT,TT,D,S}) where {DT,TT,D,S} = D +@inline Base.ndims(::IntegratorFLRK{DT,TT,D,S}) where {DT,TT,D,S} = D function initialize!(int::IntegratorFLRK, sol::AtomicSolutionODE) diff --git a/src/integrators/rk/integrators_iprk.jl b/src/integrators/rk/integrators_iprk.jl index 755c8437b..9b7e7cf1d 100644 --- a/src/integrators/rk/integrators_iprk.jl +++ b/src/integrators/rk/integrators_iprk.jl @@ -162,7 +162,7 @@ struct IntegratorIPRK{DT, TT, D, S, PT <: ParametersIPRK{DT,TT}, end -@inline Base.ndims(int::IntegratorIPRK{DT,TT,D,S}) where {DT,TT,D,S} = D +@inline Base.ndims(::IntegratorIPRK{DT,TT,D,S}) where {DT,TT,D,S} = D function initialize!(int::IntegratorIPRK, sol::AtomicSolutionPODE) diff --git a/src/integrators/rk/integrators_pglrk.jl b/src/integrators/rk/integrators_pglrk.jl index 4c6346414..f19233c61 100644 --- a/src/integrators/rk/integrators_pglrk.jl +++ b/src/integrators/rk/integrators_pglrk.jl @@ -179,8 +179,8 @@ struct IntegratorPGLRK{DT, TT, D, S, PT <: ParametersPGLRK{DT,TT}, end -@inline Base.ndims(int::IntegratorPGLRK{DT,TT,D,S}) where {DT,TT,D,S} = D -@inline nstages(integrator::IntegratorPGLRK{DT,TT,D,S}) where {DT,TT,D,S} = S +@inline Base.ndims(::IntegratorPGLRK{DT,TT,D,S}) where {DT,TT,D,S} = D +@inline nstages(::IntegratorPGLRK{DT,TT,D,S}) where {DT,TT,D,S} = S function initialize!(int::IntegratorPGLRK, sol::AtomicSolutionODE) diff --git a/src/integrators/spark/integrators_hpark.jl b/src/integrators/spark/integrators_hpark.jl index 99282d856..726490fe1 100644 --- a/src/integrators/spark/integrators_hpark.jl +++ b/src/integrators/spark/integrators_hpark.jl @@ -79,6 +79,9 @@ struct IntegratorHPARK{DT, TT, D, S, R, PT <: ParametersHPARK{DT,TT,D,S,R}, end +CommonFunctions.nconstraints(::IntegratorHPARK{DT,TT,D}) where {DT,TT,D} = D + + function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R}, params::ParametersHPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} local tqᵢ::TT diff --git a/src/integrators/spark/integrators_hspark.jl b/src/integrators/spark/integrators_hspark.jl index 7562b89f9..156e0a482 100644 --- a/src/integrators/spark/integrators_hspark.jl +++ b/src/integrators/spark/integrators_hspark.jl @@ -77,6 +77,9 @@ struct IntegratorHSPARK{DT, TT, D, S, R, PT <: ParametersHSPARK{DT,TT,D,S,R}, end +CommonFunctions.nconstraints(::IntegratorHSPARK{DT,TT,D}) where {DT,TT,D} = D + + function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R}, params::ParametersHSPARK{DT,TT,D,S,R,P}) where {ST,DT,TT,D,S,R,P} local tqᵢ::TT diff --git a/src/integrators/spark/integrators_hspark_primary.jl b/src/integrators/spark/integrators_hspark_primary.jl index 8d96b90a8..2cf1b3ac4 100644 --- a/src/integrators/spark/integrators_hspark_primary.jl +++ b/src/integrators/spark/integrators_hspark_primary.jl @@ -79,6 +79,9 @@ struct IntegratorHSPARKprimary{DT, TT, D, S, R, PT <: ParametersHSPARKprimary{DT end +CommonFunctions.nconstraints(::IntegratorHSPARKprimary{DT,TT,D}) where {DT,TT,D} = D + + function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R}, params::ParametersHSPARKprimary{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} local tqᵢ::TT diff --git a/src/integrators/spark/integrators_hspark_secondary.jl b/src/integrators/spark/integrators_hspark_secondary.jl index c29bdcc84..559142051 100644 --- a/src/integrators/spark/integrators_hspark_secondary.jl +++ b/src/integrators/spark/integrators_hspark_secondary.jl @@ -82,6 +82,9 @@ struct IntegratorHSPARKsecondary{DT, TT, D, S, R, PT <: ParametersHSPARKsecondar end +CommonFunctions.nconstraints(::IntegratorHSPARKsecondary{DT,TT,D}) where {DT,TT,D} = D + + function initial_guess!(int::IntegratorHSPARKsecondary{DT}, sol::AtomicSolutionPDAE{DT}, cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT} for i in eachstage(int) diff --git a/src/integrators/spark/integrators_slrk.jl b/src/integrators/spark/integrators_slrk.jl index 49e235085..6accc1084 100644 --- a/src/integrators/spark/integrators_slrk.jl +++ b/src/integrators/spark/integrators_slrk.jl @@ -130,6 +130,9 @@ struct IntegratorSLRK{DT, TT, D, S, PT <: ParametersSLRK{DT,TT,D,S,S}, end +CommonFunctions.nconstraints(::IntegratorSLRK{DT,TT,D,S}) where {DT,TT,D,S} = D + + function initial_guess!(int::IntegratorSLRK{DT}, sol::AtomicSolutionPDAE{DT}, cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT} for i in 1:nstages(int) diff --git a/src/integrators/spark/integrators_vpark.jl b/src/integrators/spark/integrators_vpark.jl index 16023a3a0..6850252de 100644 --- a/src/integrators/spark/integrators_vpark.jl +++ b/src/integrators/spark/integrators_vpark.jl @@ -83,6 +83,9 @@ struct IntegratorVPARK{DT, TT, D, S, R, PT <: ParametersVPARK{DT,TT,D,S,R}, end +CommonFunctions.nconstraints(::IntegratorVPARK{DT,TT,D}) where {DT,TT,D} = D + + function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R}, params::ParametersVPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} local tpᵢ::TT diff --git a/src/integrators/spark/integrators_vspark.jl b/src/integrators/spark/integrators_vspark.jl index 3b936595d..55a00f94e 100644 --- a/src/integrators/spark/integrators_vspark.jl +++ b/src/integrators/spark/integrators_vspark.jl @@ -83,6 +83,9 @@ struct IntegratorVSPARK{DT, TT, D, S, R, PT <: ParametersVSPARK{DT,TT,D,S,R}, end +CommonFunctions.nconstraints(::IntegratorVSPARK{DT,TT,D}) where {DT,TT,D} = D + + function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R}, params::ParametersVSPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R} local tpᵢ::TT diff --git a/src/integrators/spark/integrators_vspark_primary.jl b/src/integrators/spark/integrators_vspark_primary.jl index 35b8f816f..c721beb09 100644 --- a/src/integrators/spark/integrators_vspark_primary.jl +++ b/src/integrators/spark/integrators_vspark_primary.jl @@ -83,6 +83,9 @@ struct IntegratorVSPARKprimary{DT, TT, D, S, R, PT <: ParametersVSPARKprimary{DT end +CommonFunctions.nconstraints(::IntegratorVSPARKprimary{DT,TT,D}) where {DT,TT,D} = D + + function initial_guess!(int::IntegratorVSPARKprimary{DT}, sol::AtomicSolutionPDAE{DT}, cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT} for i in eachstage(int) diff --git a/src/integrators/spark/integrators_vspark_secondary.jl b/src/integrators/spark/integrators_vspark_secondary.jl index d51159e00..d6315225b 100644 --- a/src/integrators/spark/integrators_vspark_secondary.jl +++ b/src/integrators/spark/integrators_vspark_secondary.jl @@ -134,6 +134,9 @@ struct IntegratorVSPARKsecondary{DT, TT, D, S, R, end +CommonFunctions.nconstraints(::IntegratorVSPARKsecondary{DT,TT,D}) where {DT,TT,D} = D + + function initial_guess!(int::IntegratorVSPARKsecondary{DT}, sol::AtomicSolutionPDAE{DT}, cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT} for i in 1:pstages(int) diff --git a/src/integrators/splitting/integrators_composition.jl b/src/integrators/splitting/integrators_composition.jl index 0b4de24b2..58169f231 100644 --- a/src/integrators/splitting/integrators_composition.jl +++ b/src/integrators/splitting/integrators_composition.jl @@ -48,6 +48,7 @@ function IntegratorComposition(equation::SODE{DT}, tableau::AbstractTableauSplit end +@inline Base.ndims(::IntegratorComposition{DT,TT,D}) where {DT,TT,D} = D timestep(int::IntegratorComposition) = int.Δt diff --git a/src/integrators/splitting/integrators_exact_ode.jl b/src/integrators/splitting/integrators_exact_ode.jl index c9524a23d..60745b12f 100644 --- a/src/integrators/splitting/integrators_exact_ode.jl +++ b/src/integrators/splitting/integrators_exact_ode.jl @@ -10,6 +10,7 @@ struct IntegratorExactODE{DT, TT, D, QT <: Function} <: DeterministicIntegrator{ end +@inline Base.ndims(::IntegratorExactODE{DT,TT,D}) where {DT,TT,D} = D timestep(int::IntegratorExactODE) = int.Δt diff --git a/src/integrators/splitting/integrators_splitting.jl b/src/integrators/splitting/integrators_splitting.jl index bdcbc9673..bbb9694e6 100644 --- a/src/integrators/splitting/integrators_splitting.jl +++ b/src/integrators/splitting/integrators_splitting.jl @@ -21,6 +21,7 @@ function IntegratorSplitting(equation::SODE{DT,TT}, tableau::ST, Δt::TT) where end +@inline Base.ndims(::IntegratorSplitting{DT,TT,D}) where {DT,TT,D} = D timestep(int::IntegratorSplitting) = int.Δt From 0566527cc9dae4a3beebebdce1473d3ecff055bc Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:27:47 +0200 Subject: [PATCH 06/30] Add support for storing internal state variables of an integrator and solver status in atomic solution. (+1 squashed commit) Squashed commits: [0e5695a] Add internal state to atomic solutions for SDEs. --- src/integrators/abstract_integrator.jl | 15 +++++++++++++++ src/solutions/atomic_solution_dae.jl | 11 +++++++---- src/solutions/atomic_solution_ode.jl | 11 +++++++---- src/solutions/atomic_solution_pdae.jl | 11 +++++++---- src/solutions/atomic_solution_pode.jl | 11 +++++++---- src/solutions/atomic_solution_psde.jl | 10 ++++++---- src/solutions/atomic_solution_sde.jl | 10 ++++++---- 7 files changed, 55 insertions(+), 24 deletions(-) diff --git a/src/integrators/abstract_integrator.jl b/src/integrators/abstract_integrator.jl index 60a52a039..e1b58bb6c 100644 --- a/src/integrators/abstract_integrator.jl +++ b/src/integrators/abstract_integrator.jl @@ -27,6 +27,21 @@ CommonFunctions.nconstraints(integrator::Integrator) = error("nconstraints() not eachdim(integrator::Integrator) = 1:ndims(integrator) +get_internal_variables(::Integrator) = NamedTuple() + + +Solutions.AtomicSolution(integrator::ODEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionODE(DT, TT, ndims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::PODEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPODE(DT, TT, ndims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::DAEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::PDAEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator)) + abstract type Parameters{DT,TT} end diff --git a/src/solutions/atomic_solution_dae.jl b/src/solutions/atomic_solution_dae.jl index 6255d36a3..8d7c7b9db 100644 --- a/src/solutions/atomic_solution_dae.jl +++ b/src/solutions/atomic_solution_dae.jl @@ -15,7 +15,7 @@ Atomic solution for an DAE. * `u`: projective vector field of q * `u̅`: projective vector field of q̅ """ -mutable struct AtomicSolutionDAE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionDAE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -32,15 +32,18 @@ mutable struct AtomicSolutionDAE{DT,TT} <: AtomicSolution{DT,TT} u::Vector{DT} u̅::Vector{DT} - function AtomicSolutionDAE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionDAE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nm), zeros(DT, nm), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), + internal) end end -AtomicSolutionDAE(DT, TT, nd, nm) = AtomicSolutionDAE{DT, TT}(nd, nm) +AtomicSolutionDAE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionDAE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionDAE, sol) t, q, λ = sol diff --git a/src/solutions/atomic_solution_ode.jl b/src/solutions/atomic_solution_ode.jl index 476e11df6..5292bc1db 100644 --- a/src/solutions/atomic_solution_ode.jl +++ b/src/solutions/atomic_solution_ode.jl @@ -11,7 +11,7 @@ Atomic solution for an ODE. * `v`: vector field of q * `v̅`: vector field of q̅ """ -mutable struct AtomicSolutionODE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionODE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -22,13 +22,16 @@ mutable struct AtomicSolutionODE{DT,TT} <: AtomicSolution{DT,TT} v::Vector{DT} v̅::Vector{DT} - function AtomicSolutionODE{DT, TT}(nd) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionODE{DT, TT, IT}(nd, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), internal) end end -AtomicSolutionODE(DT, TT, nd) = AtomicSolutionODE{DT, TT}(nd) +AtomicSolutionODE(DT, TT, nd, internal::IT=NamedTuple()) where {IT} = AtomicSolutionODE{DT, TT, IT}(nd, internal) + function set_solution!(asol::AtomicSolutionODE, sol) t, q = sol diff --git a/src/solutions/atomic_solution_pdae.jl b/src/solutions/atomic_solution_pdae.jl index 3488a76b7..0d4e2ab5f 100644 --- a/src/solutions/atomic_solution_pdae.jl +++ b/src/solutions/atomic_solution_pdae.jl @@ -22,7 +22,7 @@ Atomic solution for an PDAE. * `g`: projective vector field of p * `g̅`: projective vector field of p̅ """ -mutable struct AtomicSolutionPDAE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionPDAE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -47,16 +47,19 @@ mutable struct AtomicSolutionPDAE{DT,TT} <: AtomicSolution{DT,TT} g::Vector{DT} g̅::Vector{DT} - function AtomicSolutionPDAE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionPDAE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nm), zeros(DT, nm), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), + internal) end end -AtomicSolutionPDAE(DT, TT, nd, nm) = AtomicSolutionPDAE{DT, TT}(nd, nm) +AtomicSolutionPDAE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionPDAE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionPDAE, sol) t, q, p, λ = sol diff --git a/src/solutions/atomic_solution_pode.jl b/src/solutions/atomic_solution_pode.jl index 53aee9cee..31b56eab9 100644 --- a/src/solutions/atomic_solution_pode.jl +++ b/src/solutions/atomic_solution_pode.jl @@ -16,7 +16,7 @@ Atomic solution for an PODE. * `f`: vector field of p * `f̅`: vector field of p̅ """ -mutable struct AtomicSolutionPODE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionPODE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -33,14 +33,17 @@ mutable struct AtomicSolutionPODE{DT,TT} <: AtomicSolution{DT,TT} f::Vector{DT} f̅::Vector{DT} - function AtomicSolutionPODE{DT, TT}(nd) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionPODE{DT, TT, IT}(nd, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), + internal) end end -AtomicSolutionPODE(DT, TT, nd) = AtomicSolutionPODE{DT, TT}(nd) +AtomicSolutionPODE(DT, TT, nd, internal::IT=NamedTuple()) where {IT} = AtomicSolutionPODE{DT, TT, IT}(nd, internal) function set_solution!(asol::AtomicSolutionPODE, sol) t, q, p = sol diff --git a/src/solutions/atomic_solution_psde.jl b/src/solutions/atomic_solution_psde.jl index 19d27d7d7..3c9360980 100644 --- a/src/solutions/atomic_solution_psde.jl +++ b/src/solutions/atomic_solution_psde.jl @@ -15,7 +15,7 @@ Atomic solution for an SDE. * `ΔZ`: Wiener process driving the stochastic process q * `K`: integer parameter defining the truncation of the increments of the Wiener process (for strong solutions) """ -mutable struct AtomicSolutionPSDE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionPSDE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -31,14 +31,16 @@ mutable struct AtomicSolutionPSDE{DT,TT} <: AtomicSolution{DT,TT} ΔZ::Vector{DT} K::Int - function AtomicSolutionPSDE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionPSDE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nm), zeros(DT, nm)) + zeros(DT, nm), zeros(DT, nm), 0, internal) end end -AtomicSolutionPSDE(DT, TT, nd, nm) = AtomicSolutionPSDE{DT, TT}(nd, nm) +AtomicSolutionPSDE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionPSDE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionPSDE, sol) t, q, p = sol diff --git a/src/solutions/atomic_solution_sde.jl b/src/solutions/atomic_solution_sde.jl index b3292f704..1bab63949 100644 --- a/src/solutions/atomic_solution_sde.jl +++ b/src/solutions/atomic_solution_sde.jl @@ -12,7 +12,7 @@ Atomic solution for an SDE. * `ΔZ`: Wiener process driving the stochastic process q * `K`: integer parameter defining the truncation of the increments of the Wiener process (for strong solutions) """ -mutable struct AtomicSolutionSDE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionSDE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -24,13 +24,15 @@ mutable struct AtomicSolutionSDE{DT,TT} <: AtomicSolution{DT,TT} ΔZ::Vector{DT} K::Int - function AtomicSolutionSDE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionSDE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nm), zeros(DT, nm)) + zeros(DT, nm), zeros(DT, nm), 0, internal) end end -AtomicSolutionSDE(DT, TT, nd, nm) = AtomicSolutionSDE{DT, TT}(nd, nm) +AtomicSolutionSDE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionSDE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionSDE, sol) t, q = sol From b009c1ac6ead753d53a6c714630b3c2598d185db Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:27:00 +0200 Subject: [PATCH 07/30] Revise integrator type hierarchy. --- src/Integrators.jl | 11 +++++-- src/integrators/SPARK.jl | 5 ++-- src/integrators/Stochastic.jl | 3 +- src/integrators/VPRK.jl | 5 ++-- src/integrators/abstract_integrator.jl | 19 +++++++++--- src/integrators/cgvi/integrators_cgvi.jl | 2 +- src/integrators/dgvi/integrators_dgvi.jl | 2 +- .../dgvi/integrators_dgvi_experimental.jl | 2 +- .../dgvi/integrators_dgvi_path_integral.jl | 2 +- .../dgvi/integrators_dgvi_projection_final.jl | 2 +- .../integrators_dgvi_projection_initial.jl | 2 +- src/integrators/integrators.jl | 5 +--- src/integrators/rk/abstract_integrator_rk.jl | 18 ++++++----- src/integrators/rk/integrators_flrk.jl | 7 ++--- src/integrators/rk/integrators_pglrk.jl | 2 +- .../spark/abstract_integrator_spark.jl | 3 +- .../splitting/integrators_composition.jl | 2 +- .../splitting/integrators_exact_ode.jl | 2 +- .../splitting/integrators_splitting.jl | 2 +- .../stochastic/abstract_stochastic_rk.jl | 30 ++++++++++--------- .../stochastic/integrators_siprk.jl | 2 +- .../stochastic/integrators_sisprk.jl | 2 +- .../vprk/integrators_vprk_abstract.jl | 2 +- test/integrators/special_integrators_tests.jl | 4 --- 24 files changed, 74 insertions(+), 62 deletions(-) diff --git a/src/Integrators.jl b/src/Integrators.jl index a5621ede8..c3f7039d9 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -33,8 +33,15 @@ module Integrators include("integrators/abstract_tableau.jl") - export Integrator, DeterministicIntegrator, StochasticIntegrator, IntegratorCache, - IntegratorConstructor + export Integrator, DeterministicIntegrator, StochasticIntegrator + export ODEIntegrator, DAEIntegrator, SDEIntegrator, + PODEIntegrator, PDAEIntegrator, PSDEIntegrator, + IODEIntegrator, IDAEIntegrator, + HODEIntegrator, HDAEIntegrator, + VODEIntegrator, VDAEIntegrator, + SPSDEIntegrator + + export IntegratorCache, IntegratorConstructor export integrate, integrate!, integrate_step!, equation, timestep export NonlinearFunctionParameters, function_stages! diff --git a/src/integrators/SPARK.jl b/src/integrators/SPARK.jl index a8e4ff6f3..4c332bb75 100644 --- a/src/integrators/SPARK.jl +++ b/src/integrators/SPARK.jl @@ -13,9 +13,8 @@ module SPARK import ..Integrators - import ..Integrators: DeterministicIntegrator, IntegratorPRK, Parameters, - IDAEIntegratorCache, InitialGuessIODE, InitialGuessPODE - import ..Integrators: IntegratorCache, CacheDict, CacheType + import ..Integrators: PDAEIntegrator, InitialGuessIODE, InitialGuessPODE, Parameters + import ..Integrators: IDAEIntegratorCache, IntegratorCache, CacheDict, CacheType import ..Integrators: AbstractTableau, AbstractTableauERK, AbstractTableauIRK, AbstractCoefficients, CoefficientsRK, @CoefficientsRK, @HeaderCoefficientsRK diff --git a/src/integrators/Stochastic.jl b/src/integrators/Stochastic.jl index 6dd7207ca..c598840ae 100644 --- a/src/integrators/Stochastic.jl +++ b/src/integrators/Stochastic.jl @@ -9,12 +9,13 @@ module Stochastic using ..Utils import ..Integrators + import ..Integrators: nstages, noisedims import ..Equations: SDE, PSDE, SPSDE, get_function_tuple import ..Solutions: AtomicSolutionSDE, AtomicSolutionPSDE, SolutionVector import ..Solutions: update! - import ..Integrators: StochasticIntegrator, Parameters + import ..Integrators: StochasticIntegrator, SDEIntegrator, PSDEIntegrator, Parameters import ..Integrators: SDEIntegratorCache, PSDEIntegratorCache, IntegratorCache, CacheDict, CacheType import ..Integrators: CoefficientsRK, AbstractTableauERK, AbstractTableauIRK diff --git a/src/integrators/VPRK.jl b/src/integrators/VPRK.jl index 10133cafb..2cd254381 100644 --- a/src/integrators/VPRK.jl +++ b/src/integrators/VPRK.jl @@ -14,9 +14,8 @@ module VPRK import ..Integrators - import ..Integrators: DeterministicIntegrator, IntegratorPRK, Parameters, - IODEIntegratorCache, InitialGuessIODE - import ..Integrators: IntegratorCache, CacheDict, CacheType + import ..Integrators: IODEIntegrator, IODEIntegratorCache, InitialGuessIODE, IntegratorPRK + import ..Integrators: IntegratorCache, CacheDict, CacheType, Parameters import ..Integrators: AbstractTableauPRK, AbstractCoefficients, CoefficientsRK, CoefficientsPGLRK, @CoefficientsRK, @HeaderTableau, @HeaderCoefficientsRK, diff --git a/src/integrators/abstract_integrator.jl b/src/integrators/abstract_integrator.jl index e1b58bb6c..7874c37fe 100644 --- a/src/integrators/abstract_integrator.jl +++ b/src/integrators/abstract_integrator.jl @@ -6,11 +6,11 @@ abstract type StochasticIntegrator{dType, tType} <: Integrator{dType, tType} end abstract type ODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end abstract type DAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end -abstract type IODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end -abstract type IDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end abstract type PODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end abstract type PDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end +abstract type IODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end +abstract type IDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end abstract type HODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end abstract type HDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end abstract type VODEIntegrator{dType, tType} <: IODEIntegrator{dType, tType} end @@ -24,6 +24,8 @@ equation(integrator::Integrator) = error("equation() not implemented for ", type timestep(integrator::Integrator) = error("timestep() not implemented for ", typeof(integrator)) Base.ndims(integrator::Integrator) = error("ndims() not implemented for ", typeof(integrator)) CommonFunctions.nconstraints(integrator::Integrator) = error("nconstraints() not implemented for ", typeof(integrator)) +noisedims(integrator::Integrator) = error("noisedims() not implemented for ", typeof(integrator)) +nstages(integrator::Integrator) = error("nstages() not implemented for ", typeof(integrator)) eachdim(integrator::Integrator) = 1:ndims(integrator) @@ -42,11 +44,20 @@ Solutions.AtomicSolution(integrator::DAEIntegrator{DT,TT}) where {DT,TT} = Solutions.AtomicSolution(integrator::PDAEIntegrator{DT,TT}) where {DT,TT} = AtomicSolutionPDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator)) +Solutions.AtomicSolution(integrator::SDEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::PSDEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::SPSDEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator)) + abstract type Parameters{DT,TT} end -function_stages!(x::Vector{DT}, b::Vector{DT}, params::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT) -solution_stages!(x::Vector{DT}, y::Vector{DT}, params::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT) +function_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT) +solution_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT) initialize!(::Integrator, ::AtomicSolution) = nothing diff --git a/src/integrators/cgvi/integrators_cgvi.jl b/src/integrators/cgvi/integrators_cgvi.jl index 93a6649c3..04a11afc3 100644 --- a/src/integrators/cgvi/integrators_cgvi.jl +++ b/src/integrators/cgvi/integrators_cgvi.jl @@ -127,7 +127,7 @@ struct IntegratorCGVI{DT, TT, D, S, R, BT <: Basis, PT <: ParametersCGVI{DT,TT,D,S,R}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessIODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessIODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi.jl b/src/integrators/dgvi/integrators_dgvi.jl index 85f82b008..5e08f86b5 100644 --- a/src/integrators/dgvi/integrators_dgvi.jl +++ b/src/integrators/dgvi/integrators_dgvi.jl @@ -380,7 +380,7 @@ struct IntegratorDGVI{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVI{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi_experimental.jl b/src/integrators/dgvi/integrators_dgvi_experimental.jl index cf5813cb9..fbe5a4a54 100644 --- a/src/integrators/dgvi/integrators_dgvi_experimental.jl +++ b/src/integrators/dgvi/integrators_dgvi_experimental.jl @@ -124,7 +124,7 @@ struct IntegratorDGVIEXP{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVIEXP{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi_path_integral.jl b/src/integrators/dgvi/integrators_dgvi_path_integral.jl index 7659c3c36..a13e87151 100644 --- a/src/integrators/dgvi/integrators_dgvi_path_integral.jl +++ b/src/integrators/dgvi/integrators_dgvi_path_integral.jl @@ -266,7 +266,7 @@ struct IntegratorDGVIPI{DT, TT, D, S, R, JT <: Discontinuity, PT <: ParametersDGVIPI{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} jump::JT diff --git a/src/integrators/dgvi/integrators_dgvi_projection_final.jl b/src/integrators/dgvi/integrators_dgvi_projection_final.jl index fdf665301..923844896 100644 --- a/src/integrators/dgvi/integrators_dgvi_projection_final.jl +++ b/src/integrators/dgvi/integrators_dgvi_projection_final.jl @@ -124,7 +124,7 @@ struct IntegratorDGVIP1{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVIP1{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi_projection_initial.jl b/src/integrators/dgvi/integrators_dgvi_projection_initial.jl index 7bea8d150..29dc35bd2 100644 --- a/src/integrators/dgvi/integrators_dgvi_projection_initial.jl +++ b/src/integrators/dgvi/integrators_dgvi_projection_initial.jl @@ -127,7 +127,7 @@ struct IntegratorDGVIP0{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVIP0{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/integrators.jl b/src/integrators/integrators.jl index 08e8a6616..623f54e73 100644 --- a/src/integrators/integrators.jl +++ b/src/integrators/integrators.jl @@ -213,7 +213,7 @@ function integrate!(int::Integrator{DT,TT}, sol::Solution{DT,TT}, m1::Int, m2::I @assert n2 ≥ n1 @assert n2 ≤ ntime(sol) - asol = AtomicSolution(sol) + asol = AtomicSolution(int) # loop over initial conditions showing progress bar for m in m1:m2 @@ -264,6 +264,3 @@ function integrate!(int::Integrator{DT,TT}, sol::Solution{DT,TT}, asol::AtomicSo # copy solution from cache to solution set_solution!(sol, asol, n, m) end - - -# TODO Add solver status information to all integrators (if requested). diff --git a/src/integrators/rk/abstract_integrator_rk.jl b/src/integrators/rk/abstract_integrator_rk.jl index 575e41246..919510f72 100644 --- a/src/integrators/rk/abstract_integrator_rk.jl +++ b/src/integrators/rk/abstract_integrator_rk.jl @@ -1,11 +1,13 @@ -abstract type IntegratorRK{dType, tType} <: DeterministicIntegrator{dType, tType} end -abstract type IntegratorPRK{dType, tType} <: IntegratorRK{dType, tType} end +abstract type IntegratorRK{dType, tType} <: ODEIntegrator{dType, tType} end +abstract type IntegratorPRK{dType, tType} <: PODEIntegrator{dType, tType} end -@inline equation(integrator::IntegratorRK, i::Symbol) = integrator.params.equs[i] -@inline equations(integrator::IntegratorRK) = integrator.params.equs -@inline timestep(integrator::IntegratorRK) = integrator.params.Δt -@inline tableau(integrator::IntegratorRK) = integrator.params.tab +AbstractIntegratorRK = Union{IntegratorRK,IntegratorPRK} -@inline nstages(integrator::IntegratorRK) = nstages(tableau(integrator)) -@inline eachstage(integrator::IntegratorRK) = 1:nstages(integrator) +@inline equation(integrator::AbstractIntegratorRK, i::Symbol) = integrator.params.equs[i] +@inline equations(integrator::AbstractIntegratorRK) = integrator.params.equs +@inline timestep(integrator::AbstractIntegratorRK) = integrator.params.Δt +@inline tableau(integrator::AbstractIntegratorRK) = integrator.params.tab + +@inline nstages(integrator::AbstractIntegratorRK) = nstages(tableau(integrator)) +@inline eachstage(integrator::AbstractIntegratorRK) = 1:nstages(integrator) diff --git a/src/integrators/rk/integrators_flrk.jl b/src/integrators/rk/integrators_flrk.jl index 47aa272a5..38d61e5f2 100644 --- a/src/integrators/rk/integrators_flrk.jl +++ b/src/integrators/rk/integrators_flrk.jl @@ -45,7 +45,7 @@ end "Formal Lagrangian Runge-Kutta integrator." struct IntegratorFLRK{DT, TT, D, S, PT <: ParametersFLRK{DT,TT}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: IntegratorRK{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IntegratorPRK{DT,TT} params::PT solver::ST iguess::IT @@ -97,9 +97,6 @@ struct IntegratorFLRK{DT, TT, D, S, PT <: ParametersFLRK{DT,TT}, IntegratorFLRK{DT, equation.d}(get_function_tuple(equation), tableau, Δt; kwargs...) end - function IntegratorFLRK(equation::ODE{DT,TT}, tableau::TableauFIRK{TT}, Δt::TT; kwargs...) where {DT,TT} - IntegratorFLRK{DT, ndims(equation)}(get_function_tuple(equation), tableau, Δt; kwargs...) - end end @@ -210,7 +207,7 @@ function integrate_step!(int::IntegratorFLRK, sol::AtomicSolutionPODE) integrate_diag_flrk!(int, sol) end -function integrate_step_flrk!(int::IntegratorFLRK{DT,TT}, sol::Union{AtomicSolutionODE{DT,TT},AtomicSolutionPODE{DT,TT}}, +function integrate_step_flrk!(int::IntegratorFLRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, cache::IntegratorCacheFLRK{DT}=int.caches[DT]) where {DT,TT} # update nonlinear solver parameters from atomic solution update_params!(int, sol) diff --git a/src/integrators/rk/integrators_pglrk.jl b/src/integrators/rk/integrators_pglrk.jl index f19233c61..69da9c55c 100644 --- a/src/integrators/rk/integrators_pglrk.jl +++ b/src/integrators/rk/integrators_pglrk.jl @@ -139,7 +139,7 @@ Projected Gauss-Legendre Runge-Kutta integrator. """ struct IntegratorPGLRK{DT, TT, D, S, PT <: ParametersPGLRK{DT,TT}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: IntegratorPRK{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IntegratorRK{DT,TT} params::PT solver::ST iguess::IT diff --git a/src/integrators/spark/abstract_integrator_spark.jl b/src/integrators/spark/abstract_integrator_spark.jl index 1a896aaa2..0f54e86fa 100644 --- a/src/integrators/spark/abstract_integrator_spark.jl +++ b/src/integrators/spark/abstract_integrator_spark.jl @@ -1,5 +1,5 @@ -abstract type AbstractIntegratorSPARK{dType, tType, D, S, R} <: IntegratorPRK{dType, tType} end +abstract type AbstractIntegratorSPARK{dType, tType, D, S, R} <: PDAEIntegrator{dType, tType} end abstract type AbstractIntegratorHSPARK{dType, tType, D, S, R} <: AbstractIntegratorSPARK{dType, tType, D, S, R} end abstract type AbstractIntegratorVSPARK{dType, tType, D, S, R} <: AbstractIntegratorSPARK{dType, tType, D, S, R} end @@ -7,3 +7,4 @@ abstract type AbstractIntegratorVSPARK{dType, tType, D, S, R} <: AbstractIntegra @inline nstages(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT,TT,D,S,R} = S @inline pstages(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT,TT,D,S,R} = R @inline Base.ndims(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT,TT,D,S,R} = D +@inline eachstage(integrator::AbstractIntegratorSPARK) = 1:nstages(integrator) diff --git a/src/integrators/splitting/integrators_composition.jl b/src/integrators/splitting/integrators_composition.jl index 58169f231..ff820002c 100644 --- a/src/integrators/splitting/integrators_composition.jl +++ b/src/integrators/splitting/integrators_composition.jl @@ -1,6 +1,6 @@ "Composition integrator." -struct IntegratorComposition{DT, TT, D, S, IT <: Tuple} <: DeterministicIntegrator{DT,TT} +struct IntegratorComposition{DT, TT, D, S, IT <: Tuple} <: ODEIntegrator{DT,TT} ints::IT Δt::TT q̅::Vector{DT} diff --git a/src/integrators/splitting/integrators_exact_ode.jl b/src/integrators/splitting/integrators_exact_ode.jl index 60745b12f..3a033edd3 100644 --- a/src/integrators/splitting/integrators_exact_ode.jl +++ b/src/integrators/splitting/integrators_exact_ode.jl @@ -1,6 +1,6 @@ "Exact solution of an ODE." -struct IntegratorExactODE{DT, TT, D, QT <: Function} <: DeterministicIntegrator{DT,TT} +struct IntegratorExactODE{DT, TT, D, QT <: Function} <: ODEIntegrator{DT,TT} q::QT Δt::TT diff --git a/src/integrators/splitting/integrators_splitting.jl b/src/integrators/splitting/integrators_splitting.jl index bbb9694e6..2d0497a3c 100644 --- a/src/integrators/splitting/integrators_splitting.jl +++ b/src/integrators/splitting/integrators_splitting.jl @@ -1,6 +1,6 @@ "Splitting integrator." -struct IntegratorSplitting{DT, TT, D, S, QT <: Tuple} <: DeterministicIntegrator{DT,TT} +struct IntegratorSplitting{DT, TT, D, S, QT <: Tuple} <: ODEIntegrator{DT,TT} q::QT f::NTuple{S,Int64} c::NTuple{S,TT} diff --git a/src/integrators/stochastic/abstract_stochastic_rk.jl b/src/integrators/stochastic/abstract_stochastic_rk.jl index 58aaac7f5..d1511892c 100644 --- a/src/integrators/stochastic/abstract_stochastic_rk.jl +++ b/src/integrators/stochastic/abstract_stochastic_rk.jl @@ -1,18 +1,20 @@ -abstract type StochasticIntegratorRK{dType, tType, D, M, S} <: StochasticIntegrator{dType, tType} end -abstract type StochasticIntegratorPRK{dType, tType, D, M, S} <: StochasticIntegratorRK{dType, tType, D, M, S} end +abstract type StochasticIntegratorRK{dType, tType, D, M, S} <: SDEIntegrator{dType, tType} end +abstract type StochasticIntegratorPRK{dType, tType, D, M, S} <: PSDEIntegrator{dType, tType} end -@inline parameters(integrator::StochasticIntegratorRK) = integrator.params -@inline equation(integrator::StochasticIntegratorRK, i::Symbol) = integrator.params.equs[i] -@inline equations(integrator::StochasticIntegratorRK) = integrator.params.equs -@inline timestep(integrator::StochasticIntegratorRK) = integrator.params.Δt -@inline tableau(integrator::StochasticIntegratorRK) = integrator.params.tab +AbstractStochasticIntegratorRK{DT,TT,D,M,S} = Union{StochasticIntegratorRK{DT,TT,D,M,S}, StochasticIntegratorPRK{DT,TT,D,M,S}} -@inline nstages(integrator::StochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = S -@inline noisedims(integrator::StochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = M -@inline Base.ndims(integrator::StochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = D -@inline Base.eltype(integrator::StochasticIntegratorRK{DT}) where {DT} = DT +@inline parameters(integrator::AbstractStochasticIntegratorRK) = integrator.params +@inline equation(integrator::AbstractStochasticIntegratorRK, i::Symbol) = integrator.params.equs[i] +@inline equations(integrator::AbstractStochasticIntegratorRK) = integrator.params.equs +@inline timestep(integrator::AbstractStochasticIntegratorRK) = integrator.params.Δt +@inline tableau(integrator::AbstractStochasticIntegratorRK) = integrator.params.tab -@inline eachstage(integrator::StochasticIntegratorRK) = 1:nstages(integrator) -@inline eachnoise(integrator::StochasticIntegratorRK) = 1:noisedims(integrator) -@inline eachdim(integrator::StochasticIntegratorRK) = 1:ndims(integrator) +@inline Integrators.nstages(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = S +@inline Integrators.noisedims(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = M +@inline Base.ndims(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = D +@inline Base.eltype(::AbstractStochasticIntegratorRK{DT}) where {DT} = DT + +@inline eachstage(integrator::AbstractStochasticIntegratorRK) = 1:nstages(integrator) +@inline eachnoise(integrator::AbstractStochasticIntegratorRK) = 1:noisedims(integrator) +@inline eachdim(integrator::AbstractStochasticIntegratorRK) = 1:ndims(integrator) diff --git a/src/integrators/stochastic/integrators_siprk.jl b/src/integrators/stochastic/integrators_siprk.jl index 28ba9c629..4958d3b3c 100644 --- a/src/integrators/stochastic/integrators_siprk.jl +++ b/src/integrators/stochastic/integrators_siprk.jl @@ -152,7 +152,7 @@ end "Stochastic implicit partitioned Runge-Kutta integrator." struct IntegratorSIPRK{DT, TT, D, M, S, PT <: ParametersSIPRK{DT,TT}, - ST <: NonlinearSolver{DT}} <: StochasticIntegratorRK{DT,TT,D,M,S} + ST <: NonlinearSolver{DT}} <: StochasticIntegratorPRK{DT,TT,D,M,S} params::PT solver::ST caches::CacheDict{PT} diff --git a/src/integrators/stochastic/integrators_sisprk.jl b/src/integrators/stochastic/integrators_sisprk.jl index ce90221b6..cdaa135c6 100644 --- a/src/integrators/stochastic/integrators_sisprk.jl +++ b/src/integrators/stochastic/integrators_sisprk.jl @@ -145,7 +145,7 @@ end "Stochastic implicit partitioned Runge-Kutta integrator." struct IntegratorSISPRK{DT, TT, D, M, S, PT <: ParametersSISPRK{DT,TT}, - ST <: NonlinearSolver{DT}} <: StochasticIntegratorRK{DT,TT,D,M,S} + ST <: NonlinearSolver{DT}} <: StochasticIntegratorPRK{DT,TT,D,M,S} params::PT solver::ST caches::CacheDict{PT} diff --git a/src/integrators/vprk/integrators_vprk_abstract.jl b/src/integrators/vprk/integrators_vprk_abstract.jl index 2f339206a..50ff1edf0 100644 --- a/src/integrators/vprk/integrators_vprk_abstract.jl +++ b/src/integrators/vprk/integrators_vprk_abstract.jl @@ -1,5 +1,5 @@ -abstract type AbstractIntegratorVPRK{DT,TT,D,S} <: DeterministicIntegrator{DT,TT} end +abstract type AbstractIntegratorVPRK{DT,TT,D,S} <: IODEIntegrator{DT,TT} end abstract type AbstractIntegratorVPRKwProjection{DT,TT,D,S} <: AbstractIntegratorVPRK{DT,TT,D,S} end @inline parameters(integrator::AbstractIntegratorVPRK) = integrator.params diff --git a/test/integrators/special_integrators_tests.jl b/test/integrators/special_integrators_tests.jl index 0583b6b16..66d5663ad 100644 --- a/test/integrators/special_integrators_tests.jl +++ b/test/integrators/special_integrators_tests.jl @@ -21,10 +21,6 @@ refx = sol.q[:,end] @testset "$(rpad("Special integrators",80))" begin - flint = IntegratorFLRK(ode, getTableauGLRK(2), Δt) - flsol = integrate(ode, flint, nt) - @test rel_err(flsol.q, refx) < 4E-12 - flint = IntegratorFLRK(vode, getTableauGLRK(2), Δt) flsol = integrate(vode, flint, nt) @test rel_err(flsol.q, refx) < 4E-12 From 37a37b8a29a0ec451192045cc673f25b8039e18d Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:28:17 +0200 Subject: [PATCH 08/30] Minor fixes in VPRK integrators. --- src/integrators/vprk/integrators_vprk_common.jl | 4 ++-- src/integrators/vprk/integrators_vprk_pinternal.jl | 2 +- src/integrators/vprk/integrators_vprk_psecondary.jl | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/integrators/vprk/integrators_vprk_common.jl b/src/integrators/vprk/integrators_vprk_common.jl index 1a88db656..5195104bf 100644 --- a/src/integrators/vprk/integrators_vprk_common.jl +++ b/src/integrators/vprk/integrators_vprk_common.jl @@ -118,7 +118,7 @@ function compute_stages_q_vprk!(Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, U: y1 += params.tab.q.a[i,j] * V[j][k] y2 += params.tab.q.â[i,j] * V[j][k] end - Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + params.pparams[:R][1] * U[1][k]) + Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * params.pparams[:R][1] * U[1][k] end end end @@ -198,7 +198,7 @@ function compute_rhs_vprk!(b::Vector{ST}, P::Vector{Vector{ST}}, F::Vector{Vecto z1 += params.tab.p.a[i,j] * (F[j][k] + R[j][k]) z2 += params.tab.p.â[i,j] * (F[j][k] + R[j][k]) end - b[D*(i-1)+k] = - (P[i][k] - params.p̅[k]) + params.Δt * (z1 + z2 + params.pparams[:R][1] * G[1][k]) + b[D*(i-1)+k] = - (P[i][k] - params.p̅[k]) + params.Δt * (z1 + z2) + params.Δt * params.pparams[:R][1] * G[1][k] end end end diff --git a/src/integrators/vprk/integrators_vprk_pinternal.jl b/src/integrators/vprk/integrators_vprk_pinternal.jl index 7716ef1a8..d7fda4039 100644 --- a/src/integrators/vprk/integrators_vprk_pinternal.jl +++ b/src/integrators/vprk/integrators_vprk_pinternal.jl @@ -98,7 +98,7 @@ function compute_projection_vprk!(x::Vector{ST}, y1 += params.tab.q.a[i,j] * V[j][k] y2 += params.tab.q.â[i,j] * V[j][k] end - Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + params.pparams[:R][1] * U[1][k]) + Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * params.pparams[:R][1] * U[1][k] end end diff --git a/src/integrators/vprk/integrators_vprk_psecondary.jl b/src/integrators/vprk/integrators_vprk_psecondary.jl index b658e6b2f..c2111eb5f 100644 --- a/src/integrators/vprk/integrators_vprk_psecondary.jl +++ b/src/integrators/vprk/integrators_vprk_psecondary.jl @@ -172,7 +172,7 @@ function compute_stages_q_vprk!(q::Vector{ST}, Q::Vector{Vector{ST}}, y3 += params.tab.q.a[i,j] * Λ[j][k] y4 += params.tab.q.â[i,j] * Λ[j][k] end - Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + y3 + y4) + Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * (y3 + y4) end end @@ -185,7 +185,7 @@ function compute_stages_q_vprk!(q::Vector{ST}, Q::Vector{Vector{ST}}, y3 += params.tab.q.b[j] * Λ[j][k] y4 += params.tab.q.b̂[j] * Λ[j][k] end - q[k] = params.q̅[k] + params.Δt * (y1 + y2 + y3 + y4) + q[k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * (y3 + y4) end end @@ -238,7 +238,7 @@ function compute_rhs_vprk!(b::Vector{ST}, P::Vector{Vector{ST}}, F::Vector{Vecto z3 += params.tab.p.a[i,j] * R[j][k] z4 += params.tab.p.â[i,j] * R[j][k] end - b[D*(i-1)+k] = (P[i][k] - params.p̅[k]) - params.Δt * (z1 + z2 + z3 + z4) + b[D*(i-1)+k] = (P[i][k] - params.p̅[k]) - params.Δt * (z1 + z2) - params.Δt * (z3 + z4) end end end @@ -277,7 +277,7 @@ function compute_rhs_vprk_projection!(b::Vector{ST}, p::Vector{ST}, # z3 += params.tab.p.b[j] * R[j][k] # z4 += params.tab.p.b̂[j] * R[j][k] # end - # b[offset+D*(S-1)+k] = (p[k] - params.p̅[k]) - params.Δt * (z1 + z2 + z3 + z4) + # b[offset+D*(S-1)+k] = (p[k] - params.p̅[k]) - params.Δt * (z1 + z2) + params.Δt * (z3 + z4) # end end From b27b2643f86b363705e352f46d60fe9498925661 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:28:39 +0200 Subject: [PATCH 09/30] Disable initial guess for Runge-Kutta methods for implicit equations. --- src/integrators/rk/integrators_firk_implicit.jl | 12 ++++++------ src/integrators/rk/integrators_srk_implicit.jl | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/integrators/rk/integrators_firk_implicit.jl b/src/integrators/rk/integrators_firk_implicit.jl index 628d238bf..a2d6918a3 100644 --- a/src/integrators/rk/integrators_firk_implicit.jl +++ b/src/integrators/rk/integrators_firk_implicit.jl @@ -127,19 +127,19 @@ function initial_guess!(int::IntegratorFIRKimplicit{DT,TT}, sol::AtomicSolutionP cache::IntegratorCacheFIRKimplicit{DT}=int.caches[DT]) where {DT,TT} # compute initial guess for internal stages - for i in eachstage(int) - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) - end + # for i in eachstage(int) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) + # end for i in eachstage(int) for k in eachdim(int) - int.solver.x[ndims(int)*(i-1)+k] = cache.V[i][k] + int.solver.x[ndims(int)*(i-1)+k] = 0#cache.V[i][k] end end # compute initial guess for solution - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) for k in eachdim(int) - int.solver.x[ndims(int)*nstages(int)+k] = cache.q[k] + int.solver.x[ndims(int)*nstages(int)+k] = sol.q[k]#cache.q[k] end end diff --git a/src/integrators/rk/integrators_srk_implicit.jl b/src/integrators/rk/integrators_srk_implicit.jl index 5b099d544..3c5ca6976 100644 --- a/src/integrators/rk/integrators_srk_implicit.jl +++ b/src/integrators/rk/integrators_srk_implicit.jl @@ -127,19 +127,19 @@ function initial_guess!(int::IntegratorSRKimplicit{DT,TT}, sol::AtomicSolutionPO cache::IntegratorCacheSRKimplicit{DT}=int.caches[DT]) where {DT,TT} # compute initial guess for internal stages - for i in eachstage(int) - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) - end + # for i in eachstage(int) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) + # end for i in eachstage(int) for k in eachdim(int) - int.solver.x[ndims(int)*(i-1)+k] = cache.V[i][k] + int.solver.x[ndims(int)*(i-1)+k] = 0#cache.V[i][k] end end # compute initial guess for solution - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) for k in eachdim(int) - int.solver.x[ndims(int)*nstages(int)+k] = cache.q[k] + int.solver.x[ndims(int)*nstages(int)+k] = sol.q[k]#cache.q[k] end end From dde769bbcaf3ab3ae0422edb9fbcabf95c92f5c3 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:29:09 +0200 Subject: [PATCH 10/30] Store internal stage in atomic solution for some VPRK integrators. --- src/integrators/vprk/integrators_vprk.jl | 20 ++++++++++++++++ .../vprk/integrators_vprk_pmidpoint.jl | 23 +++++++++++++++++++ .../vprk/integrators_vprk_psymmetric.jl | 23 +++++++++++++++++++ 3 files changed, 66 insertions(+) diff --git a/src/integrators/vprk/integrators_vprk.jl b/src/integrators/vprk/integrators_vprk.jl index 04fe77227..ee9c6c7c5 100644 --- a/src/integrators/vprk/integrators_vprk.jl +++ b/src/integrators/vprk/integrators_vprk.jl @@ -62,6 +62,17 @@ end IntegratorVPRKpNone = IntegratorVPRK +function Integrators.get_internal_variables(int::IntegratorVPRK{DT,TT,D,S}) where {DT, TT, D, S} + Q = create_internal_stage_vector(DT, D, S) + P = create_internal_stage_vector(DT, D, S) + V = create_internal_stage_vector(DT, D, S) + F = create_internal_stage_vector(DT, D, S) + + solver = get_solver_status(int.solver) + + (Q=Q, P=P, V=V, F=F, solver=solver) +end + function initial_guess!(int::IntegratorVPRK{DT}, sol::AtomicSolutionPODE{DT}, cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT} @@ -123,4 +134,13 @@ function Integrators.integrate_step!(int::IntegratorVPRK{DT,TT}, sol::AtomicSolu # copy solution to initial guess update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f) + + # copy internal stage variables + sol.internal[:Q] .= cache.Q + sol.internal[:P] .= cache.P + sol.internal[:V] .= cache.V + sol.internal[:F] .= cache.F + + # copy solver status + get_solver_status!(int.solver, sol.internal[:solver]) end diff --git a/src/integrators/vprk/integrators_vprk_pmidpoint.jl b/src/integrators/vprk/integrators_vprk_pmidpoint.jl index e5ad9e7f6..2075a8875 100644 --- a/src/integrators/vprk/integrators_vprk_pmidpoint.jl +++ b/src/integrators/vprk/integrators_vprk_pmidpoint.jl @@ -45,6 +45,19 @@ struct IntegratorVPRKpMidpoint{DT, TT, D, S, end +function Integrators.get_internal_variables(int::IntegratorVPRKpMidpoint{DT,TT,D,S}) where {DT, TT, D, S} + Q = create_internal_stage_vector(DT, D, S) + P = create_internal_stage_vector(DT, D, S) + V = create_internal_stage_vector(DT, D, S) + F = create_internal_stage_vector(DT, D, S) + λ = zeros(DT,D) + + solver = get_solver_status(int.solver) + + (Q=Q, P=P, V=V, F=F, λ=λ, solver=solver) +end + + function initial_guess!(int::IntegratorVPRKpMidpoint{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT} @@ -162,4 +175,14 @@ function Integrators.integrate_step!(int::IntegratorVPRKpMidpoint{DT,TT}, sol::A # copy solution to initial guess update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f) + + # copy internal stage variables + sol.internal.Q .= cache.Q + sol.internal.P .= cache.P + sol.internal.V .= cache.V + sol.internal.F .= cache.F + sol.internal.λ .= cache.λ + + # copy solver status + get_solver_status!(int.solver, sol.internal[:solver]) end diff --git a/src/integrators/vprk/integrators_vprk_psymmetric.jl b/src/integrators/vprk/integrators_vprk_psymmetric.jl index dfdbb20b2..f7bd5d250 100644 --- a/src/integrators/vprk/integrators_vprk_psymmetric.jl +++ b/src/integrators/vprk/integrators_vprk_psymmetric.jl @@ -45,6 +45,19 @@ struct IntegratorVPRKpSymmetric{DT, TT, D, S, end +function Integrators.get_internal_variables(int::IntegratorVPRKpSymmetric{DT,TT,D,S}) where {DT, TT, D, S} + Q = create_internal_stage_vector(DT, D, S) + P = create_internal_stage_vector(DT, D, S) + V = create_internal_stage_vector(DT, D, S) + F = create_internal_stage_vector(DT, D, S) + λ = zeros(DT,D) + + solver = get_solver_status(int.solver) + + (Q=Q, P=P, V=V, F=F, λ=λ, solver=solver) +end + + function initial_guess!(int::IntegratorVPRKpSymmetric{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT} for i in eachstage(int) @@ -158,4 +171,14 @@ function Integrators.integrate_step!(int::IntegratorVPRKpSymmetric{DT,TT}, sol:: # copy solution to initial guess update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f) + + # copy internal stage variables + sol.internal.Q .= cache.Q + sol.internal.P .= cache.P + sol.internal.V .= cache.V + sol.internal.F .= cache.F + sol.internal.λ .= cache.λ + + # copy solver status + get_solver_status!(int.solver, sol.internal[:solver]) end From 0ab45a2435428d1be636c9a431381cb2f8df976d Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 1 Sep 2020 10:50:22 +0200 Subject: [PATCH 11/30] Implement status and params function for nonlinear solvers. (+2 squashed commits) Squashed commits: [c0756ba] Export nonlinear solver params and status functions. [1b92cf9] Add tests for nonlinear solver status and params functions. --- src/Solvers.jl | 4 ++-- src/solvers/nonlinear/abstract_newton_solver.jl | 2 ++ src/solvers/nonlinear/nonlinear_solvers.jl | 2 ++ test/solvers/nonlinear_solvers_tests.jl | 4 ++++ 4 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/Solvers.jl b/src/Solvers.jl index bd9afce93..fcdf3d332 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -19,8 +19,8 @@ module Solvers export computeJacobian, check_jacobian, print_jacobian export NonlinearSolver, AbstractNewtonSolver, - NLsolveNewton, - NewtonSolver, QuasiNewtonSolver, + NLsolveNewton, NewtonSolver, QuasiNewtonSolver, + params, status, residual_initial!, residual_absolute!, residual_relative!, print_solver_status, check_solver_converged, check_solver_status, solve! diff --git a/src/solvers/nonlinear/abstract_newton_solver.jl b/src/solvers/nonlinear/abstract_newton_solver.jl index 09c0972d9..3901329fc 100644 --- a/src/solvers/nonlinear/abstract_newton_solver.jl +++ b/src/solvers/nonlinear/abstract_newton_solver.jl @@ -29,6 +29,8 @@ end function computeJacobian(s::AbstractNewtonSolver) computeJacobian(s.x, s.J, s.Jparams) end +status(solver::AbstractNewtonSolver) = solver.status +params(solver::AbstractNewtonSolver) = solver.params function check_jacobian(s::AbstractNewtonSolver) println("Condition Number of Jacobian: ", cond(s.J)) diff --git a/src/solvers/nonlinear/nonlinear_solvers.jl b/src/solvers/nonlinear/nonlinear_solvers.jl index 078115cde..9cd7d233e 100644 --- a/src/solvers/nonlinear/nonlinear_solvers.jl +++ b/src/solvers/nonlinear/nonlinear_solvers.jl @@ -4,6 +4,8 @@ using Printf abstract type NonlinearSolver{T} end solve!(s::NonlinearSolver) = error("solve! not implemented for $(typeof(s))") +status(s::NonlinearSolver) = error("status not implemented for $(typeof(s))") +params(s::NonlinearSolver) = error("params not implemented for $(typeof(s))") function solve!(s::NonlinearSolver{T}, x₀::Vector{T}) where {T} setInitialConditions!(s, x₀) diff --git a/test/solvers/nonlinear_solvers_tests.jl b/test/solvers/nonlinear_solvers_tests.jl index fa9b98b53..61f9b385d 100644 --- a/test/solvers/nonlinear_solvers_tests.jl +++ b/test/solvers/nonlinear_solvers_tests.jl @@ -18,6 +18,10 @@ end for Solver in (NewtonSolver, QuasiNewtonSolver, NLsolveNewton) x = ones(T, n) nl = Solver(x, F) + + @test params(nl) == nl.params + @test status(nl) == nl.status + solve!(nl) # println(nl.status.i, ", ", nl.status.rₐ,", ", nl.status.rᵣ,", ", nl.status.rₛ) for x in nl.x From 48ef9254d42da30f20cad7edf789c916e571ca5a Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:48:49 +0200 Subject: [PATCH 12/30] Implement get_solver_status functions for nonlinear solvers returning a dictionary with number of iterations, residuals, etc. (+1 squashed commit) Squashed commits: [8f605d5] export solver status functions. --- src/Solvers.jl | 1 + src/solvers/nonlinear/nonlinear_solvers.jl | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/src/Solvers.jl b/src/Solvers.jl index fcdf3d332..fac14b471 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -23,6 +23,7 @@ module Solvers params, status, residual_initial!, residual_absolute!, residual_relative!, print_solver_status, check_solver_converged, check_solver_status, + get_solver_status, get_solver_status!, solve! include("solvers/nonlinear/nonlinear_solvers.jl") diff --git a/src/solvers/nonlinear/nonlinear_solvers.jl b/src/solvers/nonlinear/nonlinear_solvers.jl index 9cd7d233e..d22007751 100644 --- a/src/solvers/nonlinear/nonlinear_solvers.jl +++ b/src/solvers/nonlinear/nonlinear_solvers.jl @@ -113,6 +113,26 @@ function check_solver_status(status::NonlinearSolverStatus, params::NonlinearSol end end +function get_solver_status!(status::NonlinearSolverStatus, params::NonlinearSolverParameters, status_dict::Dict) + status_dict[:nls_niter] = status.i + status_dict[:nls_atol] = status.rₐ + status_dict[:nls_rtol] = status.rᵣ + status_dict[:nls_stol] = status.rₛ + status_dict[:nls_converged] = check_solver_converged(status, params) + return status_dict +end + +get_solver_status!(solver::NonlinearSolver{T}, status_dict::Dict) where {T} = + get_solver_status!(status(solver), params(solver), status_dict) + +get_solver_status(solver::NonlinearSolver{T}) where {T} = get_solver_status!(solver, + Dict(:nls_niter => 0, + :nls_atol => zero(T), + :nls_rtol => zero(T), + :nls_stol => zero(T), + :nls_converged => false) + ) + function getLinearSolver(x::AbstractVector{T}) where {T} n = length(x) From fa26538b289d0d38dbe3f47ca00e72f49ab88e93 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:45:17 +0200 Subject: [PATCH 13/30] Minor cleanup in abstract Newton solver. --- src/solvers/nonlinear/abstract_newton_solver.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/solvers/nonlinear/abstract_newton_solver.jl b/src/solvers/nonlinear/abstract_newton_solver.jl index 3901329fc..b850a3397 100644 --- a/src/solvers/nonlinear/abstract_newton_solver.jl +++ b/src/solvers/nonlinear/abstract_newton_solver.jl @@ -21,17 +21,15 @@ abstract type AbstractNewtonSolver{T} <: NonlinearSolver{T} end status::NonlinearSolverStatus{T} end - function setInitialConditions!(s::AbstractNewtonSolver{T}, x₀::Vector{T}) where {T} s.x[:] = x₀ end -function computeJacobian(s::AbstractNewtonSolver) - computeJacobian(s.x, s.J, s.Jparams) -end status(solver::AbstractNewtonSolver) = solver.status params(solver::AbstractNewtonSolver) = solver.params +computeJacobian(s::AbstractNewtonSolver) = computeJacobian(s.x, s.J, s.Jparams) + function check_jacobian(s::AbstractNewtonSolver) println("Condition Number of Jacobian: ", cond(s.J)) println("Determinant of Jacobian: ", det(s.J)) From 564a7ce1ede3d6a48d3aeacd6195800098047e0b Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 21:46:10 +0200 Subject: [PATCH 14/30] Cleanup atomic solution constructor calls. --- src/solutions/atomic_solution.jl | 12 ++++++------ src/solutions/solution_dae.jl | 2 +- src/solutions/solution_ode.jl | 2 +- src/solutions/solution_pdae.jl | 2 +- src/solutions/solution_pode.jl | 2 +- src/solutions/solution_psde.jl | 2 +- src/solutions/solution_sde.jl | 2 +- 7 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/solutions/atomic_solution.jl b/src/solutions/atomic_solution.jl index dc1d899a8..3ca36db3c 100644 --- a/src/solutions/atomic_solution.jl +++ b/src/solutions/atomic_solution.jl @@ -7,32 +7,32 @@ abstract type AtomicSolution{dType <: Number, tType <: Real} end "Create AtomicSolution for ODE." function AtomicSolution(equation::AbstractEquationODE{DT,TT}) where {DT,TT} - AtomicSolutionODE{DT,TT}(ndims(equation)) + AtomicSolutionODE(DT, TT, ndims(equation)) end "Create AtomicSolution for partitioned ODE." function AtomicSolution(equation::AbstractEquationPODE{DT,TT}) where {DT,TT} - AtomicSolutionPODE{DT,TT}(ndims(equation)) + AtomicSolutionPODE(DT, TT, ndims(equation)) end "Create AtomicSolution for DAE." function AtomicSolution(equation::AbstractEquationDAE{DT,TT}) where {DT,TT} - AtomicSolutionDAE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionDAE(DT, TT, ndims(equation), equation.m) end "Create AtomicSolution for partitioned DAE." function AtomicSolution(equation::AbstractEquationPDAE{DT,TT}) where {DT,TT} - AtomicSolutionPDAE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionPDAE(DT, TT, ndims(equation), equation.m) end "Create AtomicSolution for SDE." function AtomicSolution(equation::AbstractEquationSDE{DT,TT}) where {DT,TT} - AtomicSolutionSDE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionSDE(DT, TT, ndims(equation), equation.m) end "Create AtomicSolution for PSDE." function AtomicSolution(equation::AbstractEquationPSDE{DT,TT}) where {DT,TT} - AtomicSolutionPSDE{DT,TT}(ndims(equation), equation.m) + AtomicSolutionPSDE(DT, TT, ndims(equation), equation.m) end "Print error for AtomicSolutions of equations not implemented, yet." diff --git a/src/solutions/solution_dae.jl b/src/solutions/solution_dae.jl index ed9bb8d1a..b1515cc40 100644 --- a/src/solutions/solution_dae.jl +++ b/src/solutions/solution_dae.jl @@ -157,7 +157,7 @@ Base.:(==)(sol1::SolutionDAE{DT1,TT1,N1}, sol2::SolutionDAE{DT2,TT2,N2}) where { "Create AtomicSolution for DAE." function AtomicSolution(solution::SolutionDAE{DT,TT}) where {DT,TT} - AtomicSolutionDAE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionDAE(DT, TT, solution.nd, solution.nm) end diff --git a/src/solutions/solution_ode.jl b/src/solutions/solution_ode.jl index 513851aa1..b94556d17 100644 --- a/src/solutions/solution_ode.jl +++ b/src/solutions/solution_ode.jl @@ -142,7 +142,7 @@ Base.:(==)(sol1::SolutionODE{DT1,TT1,N1}, sol2::SolutionODE{DT2,TT2,N2}) where { "Create AtomicSolution for ODE." function AtomicSolution(solution::SolutionODE{DT,TT}) where {DT,TT} - AtomicSolutionODE{DT,TT}(solution.nd) + AtomicSolutionODE(DT, TT, solution.nd) end diff --git a/src/solutions/solution_pdae.jl b/src/solutions/solution_pdae.jl index 09a669837..e99c254cc 100644 --- a/src/solutions/solution_pdae.jl +++ b/src/solutions/solution_pdae.jl @@ -163,7 +163,7 @@ Base.:(==)(sol1::SolutionPDAE{DT1,TT1,N1}, sol2::SolutionPDAE{DT2,TT2,N2}) where "Create AtomicSolution for partitioned DAE." function AtomicSolution(solution::SolutionPDAE{DT,TT}) where {DT,TT} - AtomicSolutionPDAE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionPDAE(DT, TT, solution.nd, solution.nm) end diff --git a/src/solutions/solution_pode.jl b/src/solutions/solution_pode.jl index 0902b200a..375100239 100644 --- a/src/solutions/solution_pode.jl +++ b/src/solutions/solution_pode.jl @@ -147,7 +147,7 @@ Base.:(==)(sol1::SolutionPODE, sol2::SolutionPODE) = ( "Create AtomicSolution for partitioned ODE." function AtomicSolution(solution::SolutionPODE{DT,TT}) where {DT,TT} - AtomicSolutionPODE{DT,TT}(solution.nd) + AtomicSolutionPODE(DT, TT, solution.nd) end diff --git a/src/solutions/solution_psde.jl b/src/solutions/solution_psde.jl index e438c4b60..2b37a5946 100644 --- a/src/solutions/solution_psde.jl +++ b/src/solutions/solution_psde.jl @@ -259,7 +259,7 @@ Base.:(==)(sol1::SolutionPSDE{DT1,TT1,NQ1,NW1,C1}, sol2::SolutionPSDE{DT2,TT2,NQ "Create AtomicSolution for PSDE." function AtomicSolution(solution::SolutionPSDE{DT,TT}) where {DT,TT} - AtomicSolutionPSDE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionPSDE(DT, TT, solution.nd, solution.nm) end diff --git a/src/solutions/solution_sde.jl b/src/solutions/solution_sde.jl index 91060d0c1..5f19ba4d0 100644 --- a/src/solutions/solution_sde.jl +++ b/src/solutions/solution_sde.jl @@ -250,7 +250,7 @@ Base.:(==)(sol1::SolutionSDE{DT1,TT1,NQ1,NW1,C1}, sol2::SolutionSDE{DT2,TT2,NQ2, "Create AtomicSolution for SDE." function AtomicSolution(solution::SolutionSDE{DT,TT}) where {DT,TT} - AtomicSolutionSDE{DT,TT}(solution.nd, solution.nm) + AtomicSolutionSDE(DT, TT, solution.nd, solution.nm) end From 2cfaad8fc50248b450a070db06a4ff370c588137 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:27:47 +0200 Subject: [PATCH 15/30] Add support for storing internal state variables of an integrator and solver status in atomic solution. (+1 squashed commit) Squashed commits: [0e5695a] Add internal state to atomic solutions for SDEs. --- src/integrators/abstract_integrator.jl | 15 +++++++++++++++ src/solutions/atomic_solution_dae.jl | 11 +++++++---- src/solutions/atomic_solution_ode.jl | 11 +++++++---- src/solutions/atomic_solution_pdae.jl | 11 +++++++---- src/solutions/atomic_solution_pode.jl | 11 +++++++---- src/solutions/atomic_solution_psde.jl | 10 ++++++---- src/solutions/atomic_solution_sde.jl | 10 ++++++---- 7 files changed, 55 insertions(+), 24 deletions(-) diff --git a/src/integrators/abstract_integrator.jl b/src/integrators/abstract_integrator.jl index 60a52a039..e1b58bb6c 100644 --- a/src/integrators/abstract_integrator.jl +++ b/src/integrators/abstract_integrator.jl @@ -27,6 +27,21 @@ CommonFunctions.nconstraints(integrator::Integrator) = error("nconstraints() not eachdim(integrator::Integrator) = 1:ndims(integrator) +get_internal_variables(::Integrator) = NamedTuple() + + +Solutions.AtomicSolution(integrator::ODEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionODE(DT, TT, ndims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::PODEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPODE(DT, TT, ndims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::DAEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::PDAEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator)) + abstract type Parameters{DT,TT} end diff --git a/src/solutions/atomic_solution_dae.jl b/src/solutions/atomic_solution_dae.jl index 6255d36a3..8d7c7b9db 100644 --- a/src/solutions/atomic_solution_dae.jl +++ b/src/solutions/atomic_solution_dae.jl @@ -15,7 +15,7 @@ Atomic solution for an DAE. * `u`: projective vector field of q * `u̅`: projective vector field of q̅ """ -mutable struct AtomicSolutionDAE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionDAE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -32,15 +32,18 @@ mutable struct AtomicSolutionDAE{DT,TT} <: AtomicSolution{DT,TT} u::Vector{DT} u̅::Vector{DT} - function AtomicSolutionDAE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionDAE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nm), zeros(DT, nm), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), + internal) end end -AtomicSolutionDAE(DT, TT, nd, nm) = AtomicSolutionDAE{DT, TT}(nd, nm) +AtomicSolutionDAE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionDAE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionDAE, sol) t, q, λ = sol diff --git a/src/solutions/atomic_solution_ode.jl b/src/solutions/atomic_solution_ode.jl index 476e11df6..5292bc1db 100644 --- a/src/solutions/atomic_solution_ode.jl +++ b/src/solutions/atomic_solution_ode.jl @@ -11,7 +11,7 @@ Atomic solution for an ODE. * `v`: vector field of q * `v̅`: vector field of q̅ """ -mutable struct AtomicSolutionODE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionODE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -22,13 +22,16 @@ mutable struct AtomicSolutionODE{DT,TT} <: AtomicSolution{DT,TT} v::Vector{DT} v̅::Vector{DT} - function AtomicSolutionODE{DT, TT}(nd) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionODE{DT, TT, IT}(nd, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), internal) end end -AtomicSolutionODE(DT, TT, nd) = AtomicSolutionODE{DT, TT}(nd) +AtomicSolutionODE(DT, TT, nd, internal::IT=NamedTuple()) where {IT} = AtomicSolutionODE{DT, TT, IT}(nd, internal) + function set_solution!(asol::AtomicSolutionODE, sol) t, q = sol diff --git a/src/solutions/atomic_solution_pdae.jl b/src/solutions/atomic_solution_pdae.jl index 3488a76b7..0d4e2ab5f 100644 --- a/src/solutions/atomic_solution_pdae.jl +++ b/src/solutions/atomic_solution_pdae.jl @@ -22,7 +22,7 @@ Atomic solution for an PDAE. * `g`: projective vector field of p * `g̅`: projective vector field of p̅ """ -mutable struct AtomicSolutionPDAE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionPDAE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -47,16 +47,19 @@ mutable struct AtomicSolutionPDAE{DT,TT} <: AtomicSolution{DT,TT} g::Vector{DT} g̅::Vector{DT} - function AtomicSolutionPDAE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionPDAE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nm), zeros(DT, nm), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), + internal) end end -AtomicSolutionPDAE(DT, TT, nd, nm) = AtomicSolutionPDAE{DT, TT}(nd, nm) +AtomicSolutionPDAE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionPDAE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionPDAE, sol) t, q, p, λ = sol diff --git a/src/solutions/atomic_solution_pode.jl b/src/solutions/atomic_solution_pode.jl index 53aee9cee..31b56eab9 100644 --- a/src/solutions/atomic_solution_pode.jl +++ b/src/solutions/atomic_solution_pode.jl @@ -16,7 +16,7 @@ Atomic solution for an PODE. * `f`: vector field of p * `f̅`: vector field of p̅ """ -mutable struct AtomicSolutionPODE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionPODE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -33,14 +33,17 @@ mutable struct AtomicSolutionPODE{DT,TT} <: AtomicSolution{DT,TT} f::Vector{DT} f̅::Vector{DT} - function AtomicSolutionPODE{DT, TT}(nd) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionPODE{DT, TT, IT}(nd, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd)) + zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), + internal) end end -AtomicSolutionPODE(DT, TT, nd) = AtomicSolutionPODE{DT, TT}(nd) +AtomicSolutionPODE(DT, TT, nd, internal::IT=NamedTuple()) where {IT} = AtomicSolutionPODE{DT, TT, IT}(nd, internal) function set_solution!(asol::AtomicSolutionPODE, sol) t, q, p = sol diff --git a/src/solutions/atomic_solution_psde.jl b/src/solutions/atomic_solution_psde.jl index 19d27d7d7..3c9360980 100644 --- a/src/solutions/atomic_solution_psde.jl +++ b/src/solutions/atomic_solution_psde.jl @@ -15,7 +15,7 @@ Atomic solution for an SDE. * `ΔZ`: Wiener process driving the stochastic process q * `K`: integer parameter defining the truncation of the increments of the Wiener process (for strong solutions) """ -mutable struct AtomicSolutionPSDE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionPSDE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -31,14 +31,16 @@ mutable struct AtomicSolutionPSDE{DT,TT} <: AtomicSolution{DT,TT} ΔZ::Vector{DT} K::Int - function AtomicSolutionPSDE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionPSDE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nm), zeros(DT, nm)) + zeros(DT, nm), zeros(DT, nm), 0, internal) end end -AtomicSolutionPSDE(DT, TT, nd, nm) = AtomicSolutionPSDE{DT, TT}(nd, nm) +AtomicSolutionPSDE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionPSDE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionPSDE, sol) t, q, p = sol diff --git a/src/solutions/atomic_solution_sde.jl b/src/solutions/atomic_solution_sde.jl index b3292f704..1bab63949 100644 --- a/src/solutions/atomic_solution_sde.jl +++ b/src/solutions/atomic_solution_sde.jl @@ -12,7 +12,7 @@ Atomic solution for an SDE. * `ΔZ`: Wiener process driving the stochastic process q * `K`: integer parameter defining the truncation of the increments of the Wiener process (for strong solutions) """ -mutable struct AtomicSolutionSDE{DT,TT} <: AtomicSolution{DT,TT} +mutable struct AtomicSolutionSDE{DT,TT,IT} <: AtomicSolution{DT,TT} t::TT t̅::TT @@ -24,13 +24,15 @@ mutable struct AtomicSolutionSDE{DT,TT} <: AtomicSolution{DT,TT} ΔZ::Vector{DT} K::Int - function AtomicSolutionSDE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real} + internal::IT + + function AtomicSolutionSDE{DT, TT, IT}(nd, nm, internal::IT) where {DT <: Number, TT <: Real, IT <: NamedTuple} new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd), - zeros(DT, nm), zeros(DT, nm)) + zeros(DT, nm), zeros(DT, nm), 0, internal) end end -AtomicSolutionSDE(DT, TT, nd, nm) = AtomicSolutionSDE{DT, TT}(nd, nm) +AtomicSolutionSDE(DT, TT, nd, nm, internal::IT=NamedTuple()) where {IT} = AtomicSolutionSDE{DT, TT, IT}(nd, nm, internal) function set_solution!(asol::AtomicSolutionSDE, sol) t, q = sol From 0a9e6d67660233a2d3434db996ff16e28b9de469 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:27:00 +0200 Subject: [PATCH 16/30] Revise integrator type hierarchy. --- src/Integrators.jl | 11 +++++-- src/integrators/SPARK.jl | 5 ++-- src/integrators/Stochastic.jl | 3 +- src/integrators/VPRK.jl | 5 ++-- src/integrators/abstract_integrator.jl | 19 +++++++++--- src/integrators/cgvi/integrators_cgvi.jl | 2 +- src/integrators/dgvi/integrators_dgvi.jl | 2 +- .../dgvi/integrators_dgvi_experimental.jl | 2 +- .../dgvi/integrators_dgvi_path_integral.jl | 2 +- .../dgvi/integrators_dgvi_projection_final.jl | 2 +- .../integrators_dgvi_projection_initial.jl | 2 +- src/integrators/integrators.jl | 5 +--- src/integrators/rk/abstract_integrator_rk.jl | 18 ++++++----- src/integrators/rk/integrators_flrk.jl | 7 ++--- src/integrators/rk/integrators_pglrk.jl | 2 +- .../spark/abstract_integrator_spark.jl | 3 +- .../splitting/integrators_composition.jl | 2 +- .../splitting/integrators_exact_ode.jl | 2 +- .../splitting/integrators_splitting.jl | 2 +- .../stochastic/abstract_stochastic_rk.jl | 30 ++++++++++--------- .../stochastic/integrators_siprk.jl | 2 +- .../stochastic/integrators_sisprk.jl | 2 +- .../vprk/integrators_vprk_abstract.jl | 2 +- test/integrators/special_integrators_tests.jl | 4 --- 24 files changed, 74 insertions(+), 62 deletions(-) diff --git a/src/Integrators.jl b/src/Integrators.jl index a5621ede8..c3f7039d9 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -33,8 +33,15 @@ module Integrators include("integrators/abstract_tableau.jl") - export Integrator, DeterministicIntegrator, StochasticIntegrator, IntegratorCache, - IntegratorConstructor + export Integrator, DeterministicIntegrator, StochasticIntegrator + export ODEIntegrator, DAEIntegrator, SDEIntegrator, + PODEIntegrator, PDAEIntegrator, PSDEIntegrator, + IODEIntegrator, IDAEIntegrator, + HODEIntegrator, HDAEIntegrator, + VODEIntegrator, VDAEIntegrator, + SPSDEIntegrator + + export IntegratorCache, IntegratorConstructor export integrate, integrate!, integrate_step!, equation, timestep export NonlinearFunctionParameters, function_stages! diff --git a/src/integrators/SPARK.jl b/src/integrators/SPARK.jl index a8e4ff6f3..4c332bb75 100644 --- a/src/integrators/SPARK.jl +++ b/src/integrators/SPARK.jl @@ -13,9 +13,8 @@ module SPARK import ..Integrators - import ..Integrators: DeterministicIntegrator, IntegratorPRK, Parameters, - IDAEIntegratorCache, InitialGuessIODE, InitialGuessPODE - import ..Integrators: IntegratorCache, CacheDict, CacheType + import ..Integrators: PDAEIntegrator, InitialGuessIODE, InitialGuessPODE, Parameters + import ..Integrators: IDAEIntegratorCache, IntegratorCache, CacheDict, CacheType import ..Integrators: AbstractTableau, AbstractTableauERK, AbstractTableauIRK, AbstractCoefficients, CoefficientsRK, @CoefficientsRK, @HeaderCoefficientsRK diff --git a/src/integrators/Stochastic.jl b/src/integrators/Stochastic.jl index 6dd7207ca..c598840ae 100644 --- a/src/integrators/Stochastic.jl +++ b/src/integrators/Stochastic.jl @@ -9,12 +9,13 @@ module Stochastic using ..Utils import ..Integrators + import ..Integrators: nstages, noisedims import ..Equations: SDE, PSDE, SPSDE, get_function_tuple import ..Solutions: AtomicSolutionSDE, AtomicSolutionPSDE, SolutionVector import ..Solutions: update! - import ..Integrators: StochasticIntegrator, Parameters + import ..Integrators: StochasticIntegrator, SDEIntegrator, PSDEIntegrator, Parameters import ..Integrators: SDEIntegratorCache, PSDEIntegratorCache, IntegratorCache, CacheDict, CacheType import ..Integrators: CoefficientsRK, AbstractTableauERK, AbstractTableauIRK diff --git a/src/integrators/VPRK.jl b/src/integrators/VPRK.jl index 10133cafb..2cd254381 100644 --- a/src/integrators/VPRK.jl +++ b/src/integrators/VPRK.jl @@ -14,9 +14,8 @@ module VPRK import ..Integrators - import ..Integrators: DeterministicIntegrator, IntegratorPRK, Parameters, - IODEIntegratorCache, InitialGuessIODE - import ..Integrators: IntegratorCache, CacheDict, CacheType + import ..Integrators: IODEIntegrator, IODEIntegratorCache, InitialGuessIODE, IntegratorPRK + import ..Integrators: IntegratorCache, CacheDict, CacheType, Parameters import ..Integrators: AbstractTableauPRK, AbstractCoefficients, CoefficientsRK, CoefficientsPGLRK, @CoefficientsRK, @HeaderTableau, @HeaderCoefficientsRK, diff --git a/src/integrators/abstract_integrator.jl b/src/integrators/abstract_integrator.jl index e1b58bb6c..7874c37fe 100644 --- a/src/integrators/abstract_integrator.jl +++ b/src/integrators/abstract_integrator.jl @@ -6,11 +6,11 @@ abstract type StochasticIntegrator{dType, tType} <: Integrator{dType, tType} end abstract type ODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end abstract type DAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end -abstract type IODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end -abstract type IDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end abstract type PODEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end abstract type PDAEIntegrator{dType, tType} <: DeterministicIntegrator{dType, tType} end +abstract type IODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end +abstract type IDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end abstract type HODEIntegrator{dType, tType} <: PODEIntegrator{dType, tType} end abstract type HDAEIntegrator{dType, tType} <: PDAEIntegrator{dType, tType} end abstract type VODEIntegrator{dType, tType} <: IODEIntegrator{dType, tType} end @@ -24,6 +24,8 @@ equation(integrator::Integrator) = error("equation() not implemented for ", type timestep(integrator::Integrator) = error("timestep() not implemented for ", typeof(integrator)) Base.ndims(integrator::Integrator) = error("ndims() not implemented for ", typeof(integrator)) CommonFunctions.nconstraints(integrator::Integrator) = error("nconstraints() not implemented for ", typeof(integrator)) +noisedims(integrator::Integrator) = error("noisedims() not implemented for ", typeof(integrator)) +nstages(integrator::Integrator) = error("nstages() not implemented for ", typeof(integrator)) eachdim(integrator::Integrator) = 1:ndims(integrator) @@ -42,11 +44,20 @@ Solutions.AtomicSolution(integrator::DAEIntegrator{DT,TT}) where {DT,TT} = Solutions.AtomicSolution(integrator::PDAEIntegrator{DT,TT}) where {DT,TT} = AtomicSolutionPDAE(DT, TT, ndims(integrator), nconstraints(integrator), get_internal_variables(integrator)) +Solutions.AtomicSolution(integrator::SDEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::PSDEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator)) + +Solutions.AtomicSolution(integrator::SPSDEIntegrator{DT,TT}) where {DT,TT} = + AtomicSolutionPSDE(DT, TT, ndims(integrator), noisedims(integrator), get_internal_variables(integrator)) + abstract type Parameters{DT,TT} end -function_stages!(x::Vector{DT}, b::Vector{DT}, params::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT) -solution_stages!(x::Vector{DT}, y::Vector{DT}, params::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT) +function_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("function_stages!() not implemented for ", PT) +solution_stages!(::Vector{DT}, ::Vector{DT}, ::PT) where {DT, TT, PT <: Parameters{DT,TT}} = error("solution_stages!() not implemented for ", PT) initialize!(::Integrator, ::AtomicSolution) = nothing diff --git a/src/integrators/cgvi/integrators_cgvi.jl b/src/integrators/cgvi/integrators_cgvi.jl index 93a6649c3..04a11afc3 100644 --- a/src/integrators/cgvi/integrators_cgvi.jl +++ b/src/integrators/cgvi/integrators_cgvi.jl @@ -127,7 +127,7 @@ struct IntegratorCGVI{DT, TT, D, S, R, BT <: Basis, PT <: ParametersCGVI{DT,TT,D,S,R}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessIODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessIODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi.jl b/src/integrators/dgvi/integrators_dgvi.jl index 85f82b008..5e08f86b5 100644 --- a/src/integrators/dgvi/integrators_dgvi.jl +++ b/src/integrators/dgvi/integrators_dgvi.jl @@ -380,7 +380,7 @@ struct IntegratorDGVI{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVI{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi_experimental.jl b/src/integrators/dgvi/integrators_dgvi_experimental.jl index cf5813cb9..fbe5a4a54 100644 --- a/src/integrators/dgvi/integrators_dgvi_experimental.jl +++ b/src/integrators/dgvi/integrators_dgvi_experimental.jl @@ -124,7 +124,7 @@ struct IntegratorDGVIEXP{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVIEXP{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi_path_integral.jl b/src/integrators/dgvi/integrators_dgvi_path_integral.jl index 7659c3c36..a13e87151 100644 --- a/src/integrators/dgvi/integrators_dgvi_path_integral.jl +++ b/src/integrators/dgvi/integrators_dgvi_path_integral.jl @@ -266,7 +266,7 @@ struct IntegratorDGVIPI{DT, TT, D, S, R, JT <: Discontinuity, PT <: ParametersDGVIPI{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} jump::JT diff --git a/src/integrators/dgvi/integrators_dgvi_projection_final.jl b/src/integrators/dgvi/integrators_dgvi_projection_final.jl index fdf665301..923844896 100644 --- a/src/integrators/dgvi/integrators_dgvi_projection_final.jl +++ b/src/integrators/dgvi/integrators_dgvi_projection_final.jl @@ -124,7 +124,7 @@ struct IntegratorDGVIP1{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVIP1{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/dgvi/integrators_dgvi_projection_initial.jl b/src/integrators/dgvi/integrators_dgvi_projection_initial.jl index 7bea8d150..29dc35bd2 100644 --- a/src/integrators/dgvi/integrators_dgvi_projection_initial.jl +++ b/src/integrators/dgvi/integrators_dgvi_projection_initial.jl @@ -127,7 +127,7 @@ struct IntegratorDGVIP0{DT, TT, D, S, R, BT <: Basis, PT <: ParametersDGVIP0{DT,TT,D,S}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: DeterministicIntegrator{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IODEIntegrator{DT,TT} basis::BT quadrature::Quadrature{TT,R} diff --git a/src/integrators/integrators.jl b/src/integrators/integrators.jl index 08e8a6616..623f54e73 100644 --- a/src/integrators/integrators.jl +++ b/src/integrators/integrators.jl @@ -213,7 +213,7 @@ function integrate!(int::Integrator{DT,TT}, sol::Solution{DT,TT}, m1::Int, m2::I @assert n2 ≥ n1 @assert n2 ≤ ntime(sol) - asol = AtomicSolution(sol) + asol = AtomicSolution(int) # loop over initial conditions showing progress bar for m in m1:m2 @@ -264,6 +264,3 @@ function integrate!(int::Integrator{DT,TT}, sol::Solution{DT,TT}, asol::AtomicSo # copy solution from cache to solution set_solution!(sol, asol, n, m) end - - -# TODO Add solver status information to all integrators (if requested). diff --git a/src/integrators/rk/abstract_integrator_rk.jl b/src/integrators/rk/abstract_integrator_rk.jl index 575e41246..919510f72 100644 --- a/src/integrators/rk/abstract_integrator_rk.jl +++ b/src/integrators/rk/abstract_integrator_rk.jl @@ -1,11 +1,13 @@ -abstract type IntegratorRK{dType, tType} <: DeterministicIntegrator{dType, tType} end -abstract type IntegratorPRK{dType, tType} <: IntegratorRK{dType, tType} end +abstract type IntegratorRK{dType, tType} <: ODEIntegrator{dType, tType} end +abstract type IntegratorPRK{dType, tType} <: PODEIntegrator{dType, tType} end -@inline equation(integrator::IntegratorRK, i::Symbol) = integrator.params.equs[i] -@inline equations(integrator::IntegratorRK) = integrator.params.equs -@inline timestep(integrator::IntegratorRK) = integrator.params.Δt -@inline tableau(integrator::IntegratorRK) = integrator.params.tab +AbstractIntegratorRK = Union{IntegratorRK,IntegratorPRK} -@inline nstages(integrator::IntegratorRK) = nstages(tableau(integrator)) -@inline eachstage(integrator::IntegratorRK) = 1:nstages(integrator) +@inline equation(integrator::AbstractIntegratorRK, i::Symbol) = integrator.params.equs[i] +@inline equations(integrator::AbstractIntegratorRK) = integrator.params.equs +@inline timestep(integrator::AbstractIntegratorRK) = integrator.params.Δt +@inline tableau(integrator::AbstractIntegratorRK) = integrator.params.tab + +@inline nstages(integrator::AbstractIntegratorRK) = nstages(tableau(integrator)) +@inline eachstage(integrator::AbstractIntegratorRK) = 1:nstages(integrator) diff --git a/src/integrators/rk/integrators_flrk.jl b/src/integrators/rk/integrators_flrk.jl index 47aa272a5..38d61e5f2 100644 --- a/src/integrators/rk/integrators_flrk.jl +++ b/src/integrators/rk/integrators_flrk.jl @@ -45,7 +45,7 @@ end "Formal Lagrangian Runge-Kutta integrator." struct IntegratorFLRK{DT, TT, D, S, PT <: ParametersFLRK{DT,TT}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: IntegratorRK{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IntegratorPRK{DT,TT} params::PT solver::ST iguess::IT @@ -97,9 +97,6 @@ struct IntegratorFLRK{DT, TT, D, S, PT <: ParametersFLRK{DT,TT}, IntegratorFLRK{DT, equation.d}(get_function_tuple(equation), tableau, Δt; kwargs...) end - function IntegratorFLRK(equation::ODE{DT,TT}, tableau::TableauFIRK{TT}, Δt::TT; kwargs...) where {DT,TT} - IntegratorFLRK{DT, ndims(equation)}(get_function_tuple(equation), tableau, Δt; kwargs...) - end end @@ -210,7 +207,7 @@ function integrate_step!(int::IntegratorFLRK, sol::AtomicSolutionPODE) integrate_diag_flrk!(int, sol) end -function integrate_step_flrk!(int::IntegratorFLRK{DT,TT}, sol::Union{AtomicSolutionODE{DT,TT},AtomicSolutionPODE{DT,TT}}, +function integrate_step_flrk!(int::IntegratorFLRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, cache::IntegratorCacheFLRK{DT}=int.caches[DT]) where {DT,TT} # update nonlinear solver parameters from atomic solution update_params!(int, sol) diff --git a/src/integrators/rk/integrators_pglrk.jl b/src/integrators/rk/integrators_pglrk.jl index f19233c61..69da9c55c 100644 --- a/src/integrators/rk/integrators_pglrk.jl +++ b/src/integrators/rk/integrators_pglrk.jl @@ -139,7 +139,7 @@ Projected Gauss-Legendre Runge-Kutta integrator. """ struct IntegratorPGLRK{DT, TT, D, S, PT <: ParametersPGLRK{DT,TT}, ST <: NonlinearSolver{DT}, - IT <: InitialGuessODE{DT,TT}} <: IntegratorPRK{DT,TT} + IT <: InitialGuessODE{DT,TT}} <: IntegratorRK{DT,TT} params::PT solver::ST iguess::IT diff --git a/src/integrators/spark/abstract_integrator_spark.jl b/src/integrators/spark/abstract_integrator_spark.jl index 1a896aaa2..0f54e86fa 100644 --- a/src/integrators/spark/abstract_integrator_spark.jl +++ b/src/integrators/spark/abstract_integrator_spark.jl @@ -1,5 +1,5 @@ -abstract type AbstractIntegratorSPARK{dType, tType, D, S, R} <: IntegratorPRK{dType, tType} end +abstract type AbstractIntegratorSPARK{dType, tType, D, S, R} <: PDAEIntegrator{dType, tType} end abstract type AbstractIntegratorHSPARK{dType, tType, D, S, R} <: AbstractIntegratorSPARK{dType, tType, D, S, R} end abstract type AbstractIntegratorVSPARK{dType, tType, D, S, R} <: AbstractIntegratorSPARK{dType, tType, D, S, R} end @@ -7,3 +7,4 @@ abstract type AbstractIntegratorVSPARK{dType, tType, D, S, R} <: AbstractIntegra @inline nstages(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT,TT,D,S,R} = S @inline pstages(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT,TT,D,S,R} = R @inline Base.ndims(int::AbstractIntegratorSPARK{DT,TT,D,S,R}) where {DT,TT,D,S,R} = D +@inline eachstage(integrator::AbstractIntegratorSPARK) = 1:nstages(integrator) diff --git a/src/integrators/splitting/integrators_composition.jl b/src/integrators/splitting/integrators_composition.jl index 58169f231..ff820002c 100644 --- a/src/integrators/splitting/integrators_composition.jl +++ b/src/integrators/splitting/integrators_composition.jl @@ -1,6 +1,6 @@ "Composition integrator." -struct IntegratorComposition{DT, TT, D, S, IT <: Tuple} <: DeterministicIntegrator{DT,TT} +struct IntegratorComposition{DT, TT, D, S, IT <: Tuple} <: ODEIntegrator{DT,TT} ints::IT Δt::TT q̅::Vector{DT} diff --git a/src/integrators/splitting/integrators_exact_ode.jl b/src/integrators/splitting/integrators_exact_ode.jl index 60745b12f..3a033edd3 100644 --- a/src/integrators/splitting/integrators_exact_ode.jl +++ b/src/integrators/splitting/integrators_exact_ode.jl @@ -1,6 +1,6 @@ "Exact solution of an ODE." -struct IntegratorExactODE{DT, TT, D, QT <: Function} <: DeterministicIntegrator{DT,TT} +struct IntegratorExactODE{DT, TT, D, QT <: Function} <: ODEIntegrator{DT,TT} q::QT Δt::TT diff --git a/src/integrators/splitting/integrators_splitting.jl b/src/integrators/splitting/integrators_splitting.jl index bbb9694e6..2d0497a3c 100644 --- a/src/integrators/splitting/integrators_splitting.jl +++ b/src/integrators/splitting/integrators_splitting.jl @@ -1,6 +1,6 @@ "Splitting integrator." -struct IntegratorSplitting{DT, TT, D, S, QT <: Tuple} <: DeterministicIntegrator{DT,TT} +struct IntegratorSplitting{DT, TT, D, S, QT <: Tuple} <: ODEIntegrator{DT,TT} q::QT f::NTuple{S,Int64} c::NTuple{S,TT} diff --git a/src/integrators/stochastic/abstract_stochastic_rk.jl b/src/integrators/stochastic/abstract_stochastic_rk.jl index 58aaac7f5..d1511892c 100644 --- a/src/integrators/stochastic/abstract_stochastic_rk.jl +++ b/src/integrators/stochastic/abstract_stochastic_rk.jl @@ -1,18 +1,20 @@ -abstract type StochasticIntegratorRK{dType, tType, D, M, S} <: StochasticIntegrator{dType, tType} end -abstract type StochasticIntegratorPRK{dType, tType, D, M, S} <: StochasticIntegratorRK{dType, tType, D, M, S} end +abstract type StochasticIntegratorRK{dType, tType, D, M, S} <: SDEIntegrator{dType, tType} end +abstract type StochasticIntegratorPRK{dType, tType, D, M, S} <: PSDEIntegrator{dType, tType} end -@inline parameters(integrator::StochasticIntegratorRK) = integrator.params -@inline equation(integrator::StochasticIntegratorRK, i::Symbol) = integrator.params.equs[i] -@inline equations(integrator::StochasticIntegratorRK) = integrator.params.equs -@inline timestep(integrator::StochasticIntegratorRK) = integrator.params.Δt -@inline tableau(integrator::StochasticIntegratorRK) = integrator.params.tab +AbstractStochasticIntegratorRK{DT,TT,D,M,S} = Union{StochasticIntegratorRK{DT,TT,D,M,S}, StochasticIntegratorPRK{DT,TT,D,M,S}} -@inline nstages(integrator::StochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = S -@inline noisedims(integrator::StochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = M -@inline Base.ndims(integrator::StochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = D -@inline Base.eltype(integrator::StochasticIntegratorRK{DT}) where {DT} = DT +@inline parameters(integrator::AbstractStochasticIntegratorRK) = integrator.params +@inline equation(integrator::AbstractStochasticIntegratorRK, i::Symbol) = integrator.params.equs[i] +@inline equations(integrator::AbstractStochasticIntegratorRK) = integrator.params.equs +@inline timestep(integrator::AbstractStochasticIntegratorRK) = integrator.params.Δt +@inline tableau(integrator::AbstractStochasticIntegratorRK) = integrator.params.tab -@inline eachstage(integrator::StochasticIntegratorRK) = 1:nstages(integrator) -@inline eachnoise(integrator::StochasticIntegratorRK) = 1:noisedims(integrator) -@inline eachdim(integrator::StochasticIntegratorRK) = 1:ndims(integrator) +@inline Integrators.nstages(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = S +@inline Integrators.noisedims(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = M +@inline Base.ndims(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = D +@inline Base.eltype(::AbstractStochasticIntegratorRK{DT}) where {DT} = DT + +@inline eachstage(integrator::AbstractStochasticIntegratorRK) = 1:nstages(integrator) +@inline eachnoise(integrator::AbstractStochasticIntegratorRK) = 1:noisedims(integrator) +@inline eachdim(integrator::AbstractStochasticIntegratorRK) = 1:ndims(integrator) diff --git a/src/integrators/stochastic/integrators_siprk.jl b/src/integrators/stochastic/integrators_siprk.jl index 28ba9c629..4958d3b3c 100644 --- a/src/integrators/stochastic/integrators_siprk.jl +++ b/src/integrators/stochastic/integrators_siprk.jl @@ -152,7 +152,7 @@ end "Stochastic implicit partitioned Runge-Kutta integrator." struct IntegratorSIPRK{DT, TT, D, M, S, PT <: ParametersSIPRK{DT,TT}, - ST <: NonlinearSolver{DT}} <: StochasticIntegratorRK{DT,TT,D,M,S} + ST <: NonlinearSolver{DT}} <: StochasticIntegratorPRK{DT,TT,D,M,S} params::PT solver::ST caches::CacheDict{PT} diff --git a/src/integrators/stochastic/integrators_sisprk.jl b/src/integrators/stochastic/integrators_sisprk.jl index ce90221b6..cdaa135c6 100644 --- a/src/integrators/stochastic/integrators_sisprk.jl +++ b/src/integrators/stochastic/integrators_sisprk.jl @@ -145,7 +145,7 @@ end "Stochastic implicit partitioned Runge-Kutta integrator." struct IntegratorSISPRK{DT, TT, D, M, S, PT <: ParametersSISPRK{DT,TT}, - ST <: NonlinearSolver{DT}} <: StochasticIntegratorRK{DT,TT,D,M,S} + ST <: NonlinearSolver{DT}} <: StochasticIntegratorPRK{DT,TT,D,M,S} params::PT solver::ST caches::CacheDict{PT} diff --git a/src/integrators/vprk/integrators_vprk_abstract.jl b/src/integrators/vprk/integrators_vprk_abstract.jl index 2f339206a..50ff1edf0 100644 --- a/src/integrators/vprk/integrators_vprk_abstract.jl +++ b/src/integrators/vprk/integrators_vprk_abstract.jl @@ -1,5 +1,5 @@ -abstract type AbstractIntegratorVPRK{DT,TT,D,S} <: DeterministicIntegrator{DT,TT} end +abstract type AbstractIntegratorVPRK{DT,TT,D,S} <: IODEIntegrator{DT,TT} end abstract type AbstractIntegratorVPRKwProjection{DT,TT,D,S} <: AbstractIntegratorVPRK{DT,TT,D,S} end @inline parameters(integrator::AbstractIntegratorVPRK) = integrator.params diff --git a/test/integrators/special_integrators_tests.jl b/test/integrators/special_integrators_tests.jl index 0583b6b16..66d5663ad 100644 --- a/test/integrators/special_integrators_tests.jl +++ b/test/integrators/special_integrators_tests.jl @@ -21,10 +21,6 @@ refx = sol.q[:,end] @testset "$(rpad("Special integrators",80))" begin - flint = IntegratorFLRK(ode, getTableauGLRK(2), Δt) - flsol = integrate(ode, flint, nt) - @test rel_err(flsol.q, refx) < 4E-12 - flint = IntegratorFLRK(vode, getTableauGLRK(2), Δt) flsol = integrate(vode, flint, nt) @test rel_err(flsol.q, refx) < 4E-12 From 04d1edeeb4b2855aa684c167573d4f4bfee1738f Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:28:17 +0200 Subject: [PATCH 17/30] Minor fixes in VPRK integrators. --- src/integrators/vprk/integrators_vprk_common.jl | 4 ++-- src/integrators/vprk/integrators_vprk_pinternal.jl | 2 +- src/integrators/vprk/integrators_vprk_psecondary.jl | 8 ++++---- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/integrators/vprk/integrators_vprk_common.jl b/src/integrators/vprk/integrators_vprk_common.jl index 1a88db656..5195104bf 100644 --- a/src/integrators/vprk/integrators_vprk_common.jl +++ b/src/integrators/vprk/integrators_vprk_common.jl @@ -118,7 +118,7 @@ function compute_stages_q_vprk!(Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, U: y1 += params.tab.q.a[i,j] * V[j][k] y2 += params.tab.q.â[i,j] * V[j][k] end - Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + params.pparams[:R][1] * U[1][k]) + Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * params.pparams[:R][1] * U[1][k] end end end @@ -198,7 +198,7 @@ function compute_rhs_vprk!(b::Vector{ST}, P::Vector{Vector{ST}}, F::Vector{Vecto z1 += params.tab.p.a[i,j] * (F[j][k] + R[j][k]) z2 += params.tab.p.â[i,j] * (F[j][k] + R[j][k]) end - b[D*(i-1)+k] = - (P[i][k] - params.p̅[k]) + params.Δt * (z1 + z2 + params.pparams[:R][1] * G[1][k]) + b[D*(i-1)+k] = - (P[i][k] - params.p̅[k]) + params.Δt * (z1 + z2) + params.Δt * params.pparams[:R][1] * G[1][k] end end end diff --git a/src/integrators/vprk/integrators_vprk_pinternal.jl b/src/integrators/vprk/integrators_vprk_pinternal.jl index 7716ef1a8..d7fda4039 100644 --- a/src/integrators/vprk/integrators_vprk_pinternal.jl +++ b/src/integrators/vprk/integrators_vprk_pinternal.jl @@ -98,7 +98,7 @@ function compute_projection_vprk!(x::Vector{ST}, y1 += params.tab.q.a[i,j] * V[j][k] y2 += params.tab.q.â[i,j] * V[j][k] end - Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + params.pparams[:R][1] * U[1][k]) + Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * params.pparams[:R][1] * U[1][k] end end diff --git a/src/integrators/vprk/integrators_vprk_psecondary.jl b/src/integrators/vprk/integrators_vprk_psecondary.jl index b658e6b2f..c2111eb5f 100644 --- a/src/integrators/vprk/integrators_vprk_psecondary.jl +++ b/src/integrators/vprk/integrators_vprk_psecondary.jl @@ -172,7 +172,7 @@ function compute_stages_q_vprk!(q::Vector{ST}, Q::Vector{Vector{ST}}, y3 += params.tab.q.a[i,j] * Λ[j][k] y4 += params.tab.q.â[i,j] * Λ[j][k] end - Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2 + y3 + y4) + Q[i][k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * (y3 + y4) end end @@ -185,7 +185,7 @@ function compute_stages_q_vprk!(q::Vector{ST}, Q::Vector{Vector{ST}}, y3 += params.tab.q.b[j] * Λ[j][k] y4 += params.tab.q.b̂[j] * Λ[j][k] end - q[k] = params.q̅[k] + params.Δt * (y1 + y2 + y3 + y4) + q[k] = params.q̅[k] + params.Δt * (y1 + y2) + params.Δt * (y3 + y4) end end @@ -238,7 +238,7 @@ function compute_rhs_vprk!(b::Vector{ST}, P::Vector{Vector{ST}}, F::Vector{Vecto z3 += params.tab.p.a[i,j] * R[j][k] z4 += params.tab.p.â[i,j] * R[j][k] end - b[D*(i-1)+k] = (P[i][k] - params.p̅[k]) - params.Δt * (z1 + z2 + z3 + z4) + b[D*(i-1)+k] = (P[i][k] - params.p̅[k]) - params.Δt * (z1 + z2) - params.Δt * (z3 + z4) end end end @@ -277,7 +277,7 @@ function compute_rhs_vprk_projection!(b::Vector{ST}, p::Vector{ST}, # z3 += params.tab.p.b[j] * R[j][k] # z4 += params.tab.p.b̂[j] * R[j][k] # end - # b[offset+D*(S-1)+k] = (p[k] - params.p̅[k]) - params.Δt * (z1 + z2 + z3 + z4) + # b[offset+D*(S-1)+k] = (p[k] - params.p̅[k]) - params.Δt * (z1 + z2) + params.Δt * (z3 + z4) # end end From 697b20c7fbe7c6cb3bad7a44ed67bf7b6a64feb8 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:28:39 +0200 Subject: [PATCH 18/30] Disable initial guess for Runge-Kutta methods for implicit equations. --- src/integrators/rk/integrators_firk_implicit.jl | 12 ++++++------ src/integrators/rk/integrators_srk_implicit.jl | 12 ++++++------ 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/integrators/rk/integrators_firk_implicit.jl b/src/integrators/rk/integrators_firk_implicit.jl index 628d238bf..a2d6918a3 100644 --- a/src/integrators/rk/integrators_firk_implicit.jl +++ b/src/integrators/rk/integrators_firk_implicit.jl @@ -127,19 +127,19 @@ function initial_guess!(int::IntegratorFIRKimplicit{DT,TT}, sol::AtomicSolutionP cache::IntegratorCacheFIRKimplicit{DT}=int.caches[DT]) where {DT,TT} # compute initial guess for internal stages - for i in eachstage(int) - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) - end + # for i in eachstage(int) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) + # end for i in eachstage(int) for k in eachdim(int) - int.solver.x[ndims(int)*(i-1)+k] = cache.V[i][k] + int.solver.x[ndims(int)*(i-1)+k] = 0#cache.V[i][k] end end # compute initial guess for solution - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) for k in eachdim(int) - int.solver.x[ndims(int)*nstages(int)+k] = cache.q[k] + int.solver.x[ndims(int)*nstages(int)+k] = sol.q[k]#cache.q[k] end end diff --git a/src/integrators/rk/integrators_srk_implicit.jl b/src/integrators/rk/integrators_srk_implicit.jl index 5b099d544..3c5ca6976 100644 --- a/src/integrators/rk/integrators_srk_implicit.jl +++ b/src/integrators/rk/integrators_srk_implicit.jl @@ -127,19 +127,19 @@ function initial_guess!(int::IntegratorSRKimplicit{DT,TT}, sol::AtomicSolutionPO cache::IntegratorCacheSRKimplicit{DT}=int.caches[DT]) where {DT,TT} # compute initial guess for internal stages - for i in eachstage(int) - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) - end + # for i in eachstage(int) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).q.c[i]) + # end for i in eachstage(int) for k in eachdim(int) - int.solver.x[ndims(int)*(i-1)+k] = cache.V[i][k] + int.solver.x[ndims(int)*(i-1)+k] = 0#cache.V[i][k] end end # compute initial guess for solution - evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) + # evaluate!(int.iguess, sol.q̅, sol.v̅, sol.q, sol.v, cache.q, cache.v, one(TT)) for k in eachdim(int) - int.solver.x[ndims(int)*nstages(int)+k] = cache.q[k] + int.solver.x[ndims(int)*nstages(int)+k] = sol.q[k]#cache.q[k] end end From c89cd9bb2956b983ef25d98487c84f1db0b3f231 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Mon, 31 Aug 2020 22:29:09 +0200 Subject: [PATCH 19/30] Store internal stage in atomic solution for some VPRK integrators. --- src/integrators/vprk/integrators_vprk.jl | 20 ++++++++++++++++ .../vprk/integrators_vprk_pmidpoint.jl | 23 +++++++++++++++++++ .../vprk/integrators_vprk_psymmetric.jl | 23 +++++++++++++++++++ 3 files changed, 66 insertions(+) diff --git a/src/integrators/vprk/integrators_vprk.jl b/src/integrators/vprk/integrators_vprk.jl index 04fe77227..ee9c6c7c5 100644 --- a/src/integrators/vprk/integrators_vprk.jl +++ b/src/integrators/vprk/integrators_vprk.jl @@ -62,6 +62,17 @@ end IntegratorVPRKpNone = IntegratorVPRK +function Integrators.get_internal_variables(int::IntegratorVPRK{DT,TT,D,S}) where {DT, TT, D, S} + Q = create_internal_stage_vector(DT, D, S) + P = create_internal_stage_vector(DT, D, S) + V = create_internal_stage_vector(DT, D, S) + F = create_internal_stage_vector(DT, D, S) + + solver = get_solver_status(int.solver) + + (Q=Q, P=P, V=V, F=F, solver=solver) +end + function initial_guess!(int::IntegratorVPRK{DT}, sol::AtomicSolutionPODE{DT}, cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT} @@ -123,4 +134,13 @@ function Integrators.integrate_step!(int::IntegratorVPRK{DT,TT}, sol::AtomicSolu # copy solution to initial guess update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f) + + # copy internal stage variables + sol.internal[:Q] .= cache.Q + sol.internal[:P] .= cache.P + sol.internal[:V] .= cache.V + sol.internal[:F] .= cache.F + + # copy solver status + get_solver_status!(int.solver, sol.internal[:solver]) end diff --git a/src/integrators/vprk/integrators_vprk_pmidpoint.jl b/src/integrators/vprk/integrators_vprk_pmidpoint.jl index e5ad9e7f6..2075a8875 100644 --- a/src/integrators/vprk/integrators_vprk_pmidpoint.jl +++ b/src/integrators/vprk/integrators_vprk_pmidpoint.jl @@ -45,6 +45,19 @@ struct IntegratorVPRKpMidpoint{DT, TT, D, S, end +function Integrators.get_internal_variables(int::IntegratorVPRKpMidpoint{DT,TT,D,S}) where {DT, TT, D, S} + Q = create_internal_stage_vector(DT, D, S) + P = create_internal_stage_vector(DT, D, S) + V = create_internal_stage_vector(DT, D, S) + F = create_internal_stage_vector(DT, D, S) + λ = zeros(DT,D) + + solver = get_solver_status(int.solver) + + (Q=Q, P=P, V=V, F=F, λ=λ, solver=solver) +end + + function initial_guess!(int::IntegratorVPRKpMidpoint{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT} @@ -162,4 +175,14 @@ function Integrators.integrate_step!(int::IntegratorVPRKpMidpoint{DT,TT}, sol::A # copy solution to initial guess update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f) + + # copy internal stage variables + sol.internal.Q .= cache.Q + sol.internal.P .= cache.P + sol.internal.V .= cache.V + sol.internal.F .= cache.F + sol.internal.λ .= cache.λ + + # copy solver status + get_solver_status!(int.solver, sol.internal[:solver]) end diff --git a/src/integrators/vprk/integrators_vprk_psymmetric.jl b/src/integrators/vprk/integrators_vprk_psymmetric.jl index dfdbb20b2..f7bd5d250 100644 --- a/src/integrators/vprk/integrators_vprk_psymmetric.jl +++ b/src/integrators/vprk/integrators_vprk_psymmetric.jl @@ -45,6 +45,19 @@ struct IntegratorVPRKpSymmetric{DT, TT, D, S, end +function Integrators.get_internal_variables(int::IntegratorVPRKpSymmetric{DT,TT,D,S}) where {DT, TT, D, S} + Q = create_internal_stage_vector(DT, D, S) + P = create_internal_stage_vector(DT, D, S) + V = create_internal_stage_vector(DT, D, S) + F = create_internal_stage_vector(DT, D, S) + λ = zeros(DT,D) + + solver = get_solver_status(int.solver) + + (Q=Q, P=P, V=V, F=F, λ=λ, solver=solver) +end + + function initial_guess!(int::IntegratorVPRKpSymmetric{DT,TT}, sol::AtomicSolutionPODE{DT,TT}, cache::IntegratorCacheVPRK{DT}=int.caches[DT]) where {DT,TT} for i in eachstage(int) @@ -158,4 +171,14 @@ function Integrators.integrate_step!(int::IntegratorVPRKpSymmetric{DT,TT}, sol:: # copy solution to initial guess update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f) + + # copy internal stage variables + sol.internal.Q .= cache.Q + sol.internal.P .= cache.P + sol.internal.V .= cache.V + sol.internal.F .= cache.F + sol.internal.λ .= cache.λ + + # copy solver status + get_solver_status!(int.solver, sol.internal[:solver]) end From 6d03040979b0177485aec88357b673719a6e28b6 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 1 Sep 2020 11:05:08 +0200 Subject: [PATCH 20/30] Add nconstraints method for DAE. --- src/equations/dae.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/dae.jl b/src/equations/dae.jl index af46b99f4..9a43f96f0 100644 --- a/src/equations/dae.jl +++ b/src/equations/dae.jl @@ -132,8 +132,8 @@ function Base.similar(dae::DAE, t₀::TT, q₀::AbstractArray{DT}, λ₀::Abstra DAE(dae.v, dae.u, dae.ϕ, t₀, q₀, λ₀; h=h, parameters=parameters, periodicity=periodicity) end -@inline Base.ndims(dae::DAE) = dae.d - +@inline Base.ndims(equation::DAE) = equation.d +@inline CommonFunctions.nconstraints(equation::DAE) = equation.m @inline CommonFunctions.periodicity(equation::DAE) = equation.periodicity function get_function_tuple(equation::DAE{DT,TT,VT,UT,ϕT,HT,Nothing}) where {DT, TT, VT, UT, ϕT, HT} From fd415df93721af4f497ab288bd6238d3fcc91c11 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 1 Sep 2020 11:05:26 +0200 Subject: [PATCH 21/30] Add tests for nconstraints function of DAEs. --- test/equations/deterministic_equations_tests.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/equations/deterministic_equations_tests.jl b/test/equations/deterministic_equations_tests.jl index 20c3f3043..1d3d66b3a 100644 --- a/test/equations/deterministic_equations_tests.jl +++ b/test/equations/deterministic_equations_tests.jl @@ -233,6 +233,7 @@ end dae2 = DAE(dae_eqs..., x₀, λ₀) @test ndims(dae) == 2 + @test nconstraints(dae) == 1 @test periodicity(dae) == zero(x₀) @test get_function_tuple(dae) == NamedTuple{(:v, :u, :ϕ)}(dae_eqs) @@ -256,6 +257,7 @@ end pdae2 = PDAE(pdae_eqs..., q₀, p₀, λ₀) @test ndims(pdae) == 1 + @test nconstraints(pdae) == 1 @test periodicity(pdae) == zero(q₀) @test get_function_tuple(pdae) == NamedTuple{(:v, :f, :u, :g, :ϕ)}(pdae_eqs) @@ -280,6 +282,7 @@ end idae2 = IDAE(idae_eqs..., q₀, p₀, λ₀; v=pdae_v) @test ndims(idae) == 1 + @test nconstraints(idae) == 1 @test periodicity(idae) == zero(q₀) @test get_function_tuple(idae) == NamedTuple{(:ϑ, :f, :u, :g, :ϕ, :v)}((idae_eqs..., pdae_v)) @@ -305,6 +308,7 @@ end hdae3 = HDAE(hdae_eqs..., q₀, p₀) @test ndims(hdae) == 1 + @test nconstraints(hdae) == 1 @test periodicity(hdae) == zero(q₀) @test get_function_tuple(hdae) == NamedTuple{(:v, :f, :u, :g, :u̅, :g̅, :ϕ, :ψ, :h)}(hdae_eqs) @@ -334,6 +338,7 @@ end vdae4 = VDAE(vdae_eqs..., q₀, p₀; v=iode_v) @test ndims(vdae) == 1 + @test nconstraints(vdae) == 1 @test periodicity(vdae) == zero(q₀) @test get_function_tuple(vdae) == NamedTuple{(:ϑ, :f, :g, :g̅, :ϕ, :ψ, :v)}((vdae_eqs..., iode_v)) @@ -391,6 +396,7 @@ end spdae3 = SPDAE(spdae_eqs..., q₀, p₀) @test ndims(spdae) == 1 + @test nconstraints(spdae) == 1 @test periodicity(spdae) == zero(q₀) @test get_function_tuple(spdae) == NamedTuple{(:v, :f, :ϕ, :ψ)}(spdae_eqs) From 14219e823eddca367460317daa001f268dc492e1 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 1 Sep 2020 13:25:17 +0200 Subject: [PATCH 22/30] Add some tests for RK integrators. --- test/integrators/rk_integrators_tests.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/test/integrators/rk_integrators_tests.jl b/test/integrators/rk_integrators_tests.jl index 0aad9af60..452ec6172 100644 --- a/test/integrators/rk_integrators_tests.jl +++ b/test/integrators/rk_integrators_tests.jl @@ -77,7 +77,20 @@ end @test typeof(Integrator(ode, getTableauExplicitMidpoint(), Δt)) <: IntegratorERK @test typeof(Integrator(ode, getTableauCrouzeix(), Δt)) <: IntegratorDIRK @test typeof(Integrator(ode, getTableauImplicitMidpoint(), Δt)) <: IntegratorFIRK - @test typeof(Integrator(ode, getTableauGLRK(1), Δt)) <: IntegratorFIRK + @test typeof(Integrator(pode, TableauEPRK(:eprk_midpoint, 2, getTableauExplicitMidpoint().q), Δt)) <: IntegratorEPRK + @test typeof(Integrator(pode, TableauIPRK(:iprk_midpoint, 2, getTableauImplicitMidpoint().q), Δt)) <: IntegratorIPRK + + int_erk = Integrator(ode, getTableauExplicitMidpoint(), Δt) + int_dirk = Integrator(ode, getTableauCrouzeix(), Δt) + int_firk = Integrator(ode, getTableauImplicitMidpoint(), Δt) + int_eprk = Integrator(pode, TableauEPRK(:eprk_midpoint, 2, getTableauExplicitMidpoint().q), Δt) + int_iprk = Integrator(pode, TableauIPRK(:iprk_midpoint, 2, getTableauImplicitMidpoint().q), Δt) + + @test ndims(ode) == ndims(int_erk) + @test ndims(ode) == ndims(int_dirk) + @test ndims(ode) == ndims(int_firk) + @test ndims(pode) == ndims(int_eprk) + @test ndims(pode) == ndims(int_iprk) end From b5d23b86fab9a36c0147d0a9252a5b4dabafb744 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 1 Sep 2020 13:25:43 +0200 Subject: [PATCH 23/30] Bump version number to v0.4.1 and add release notes. --- Project.toml | 2 +- docs/src/releasenotes.md | 12 ++++++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f2d291db2..82a668a63 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GeometricIntegrators" uuid = "dcce2d33-59f6-5b8d-9047-0defad88ae06" authors = ["Michael Kraus "] -version = "0.4.0" +version = "0.4.1" [deps] DecFP = "55939f99-70c6-5e9b-8bb0-5071ed7d61fd" diff --git a/docs/src/releasenotes.md b/docs/src/releasenotes.md index 28f2fcf5e..0ed9be907 100644 --- a/docs/src/releasenotes.md +++ b/docs/src/releasenotes.md @@ -1,6 +1,18 @@ # Release Notes +## 0.4.1 + +### New Features + +* Atomic solutions can now store a NamedTuple of internal variables of the integrator, including nonlinear solver output +* Output of internal variables has been added to VPRK integrators + +### Fixes + +* Revision of integrator type hierarchy + + ## 0.4.0 ### New Integrators From 47727e54c5ff00320d0b76e4e3cb8895beeb9222 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 13:36:03 +0200 Subject: [PATCH 24/30] Cleanup in some tableaus. --- src/tableaus/tableaus_vspark_primary.jl | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/tableaus/tableaus_vspark_primary.jl b/src/tableaus/tableaus_vspark_primary.jl index 6b276084b..ed274e208 100644 --- a/src/tableaus/tableaus_vspark_primary.jl +++ b/src/tableaus/tableaus_vspark_primary.jl @@ -61,8 +61,8 @@ function getTableauVSPARKSymmetricProjection(name, q::CoefficientsRK{T}, p::Coef α_p = zeros(T, p.s, 2) α_p[:,1] .= 0.5 - a_q̃ = Array(transpose(hcat(zero(q.b), q.b))) - a_p̃ = Array(transpose(hcat(zero(p.b), p.b))) + a_q̃ = vcat(zero(q.b)', q.b') + a_p̃ = vcat(zero(p.b)', p.b') α_q̃ = [[0.0 0.0] [0.5 R∞*0.5]] @@ -79,11 +79,8 @@ function getTableauVSPARKSymmetricProjection(name, q::CoefficientsRK{T}, p::Coef c_λ = [ 0.0, 1.0] d_λ = [ 0.5, 0.5] - ω_λ = zeros(T, 1, 3) - ω_λ .= [0.5 R∞*0.5 0.0] - - δ_λ = zeros(T, 1, 2) - δ_λ .= [-1.0 +1.0] + ω_λ = T[0.5 R∞*0.5 0.0] + δ_λ = T[-1.0 +1.0] if length(d) == 0 From 63816c0a20b08992d4eae405a9bdd276e5d59a32 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 13:36:29 +0200 Subject: [PATCH 25/30] Fix in Lobatto omega matrix. --- src/tableaus/coefficients_lob.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/tableaus/coefficients_lob.jl b/src/tableaus/coefficients_lob.jl index ae7bdb1d6..a4c27060a 100644 --- a/src/tableaus/coefficients_lob.jl +++ b/src/tableaus/coefficients_lob.jl @@ -349,8 +349,13 @@ function get_lobatto_d_vector(s) end function get_lobatto_ω_matrix(s) - ω = zeros(s, s+1) - ω[1:s-1,1:s] .= getCoefficientsLobIIIA(s).a[2:s,1:s] - ω[s,s+1] = 1 + as = getCoefficientsLobIIIA(s).a[2:s,1:s] + es = zeros(s) + es[s] = 1 + + Q = vcat( hcat(as, zeros(s-1)), hcat(zeros(s)', 1) ) + L = vcat( as, es' ) + ω = inv(L) * Q + return ω end From 4cfa102251b6732f33fa84d5e9bca52fc6bc5353 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 13:36:45 +0200 Subject: [PATCH 26/30] Add convenience functions for Lobatto coefficients. --- src/tableaus/coefficients_lob.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/tableaus/coefficients_lob.jl b/src/tableaus/coefficients_lob.jl index a4c27060a..ba37235f9 100644 --- a/src/tableaus/coefficients_lob.jl +++ b/src/tableaus/coefficients_lob.jl @@ -307,6 +307,30 @@ function getCoefficientsLobIIIE(s, T=Float64) end end +function getCoefficientsLobIIIF(s, T=Float64) + if s == 2 + getCoefficientsLobIIIF2(T) + elseif s == 3 + getCoefficientsLobIIIF3(T) + elseif s == 4 + getCoefficientsLobIIIF4(T) + else + @error "Lobatto IIIF Tableau with " * string(s) * " stages not implemented." + end +end + +function getCoefficientsLobIIIG(s, T=Float64) + if s == 2 + getCoefficientsLobIIIG2(T) + elseif s == 3 + getCoefficientsLobIIIG3(T) + elseif s == 4 + getCoefficientsLobIIIG4(T) + else + @error "Lobatto IIIG Tableau with " * string(s) * " stages not implemented." + end +end + function get_lobatto_interstage_coefficients(s, σ=s+1, T=Float64) if s == 1 && σ == 2 From 2118a0e07aa0a37ad97e36db86944b09ecb4bd99 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 13:37:37 +0200 Subject: [PATCH 27/30] Export convenience functions for Lobatto coefficients. --- src/Tableaus.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Tableaus.jl b/src/Tableaus.jl index c1bee913f..e18c1076f 100644 --- a/src/Tableaus.jl +++ b/src/Tableaus.jl @@ -30,7 +30,8 @@ module Tableaus getCoefficientsLobIIIG2, getCoefficientsLobIIIG3, getCoefficientsLobIIIG4 export getCoefficientsLobIII, getCoefficientsLobIIIA, getCoefficientsLobIIIB, - getCoefficientsLobIIIC, getCoefficientsLobIIID, getCoefficientsLobIIIE + getCoefficientsLobIIIC, getCoefficientsLobIIID, getCoefficientsLobIIIE, + getCoefficientsLobIIIF, getCoefficientsLobIIIG include("tableaus/coefficients_lob.jl") From 30d853fba002a2b1cdfd591a9762a5f5d11a5bae Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 13:38:11 +0200 Subject: [PATCH 28/30] Add new VSPARK tableaus. --- src/Tableaus.jl | 14 +- src/tableaus/tableaus_vspark_primary.jl | 183 ++++++++++++++++++++++++ 2 files changed, 195 insertions(+), 2 deletions(-) diff --git a/src/Tableaus.jl b/src/Tableaus.jl index e18c1076f..83feabab1 100644 --- a/src/Tableaus.jl +++ b/src/Tableaus.jl @@ -125,15 +125,25 @@ module Tableaus include("tableaus/tableaus_vpark.jl") - export getTableauVSPARKMidpointProjection, + export getTableauVSPARKLobIIIAIIIBProjection, + getTableauVSPARKMidpointProjection, getTableauVSPARKSymmetricProjection, + getTableauVSPARKSymmetricLobProjection, getTableauVSPARKSymplecticProjection, + getTableauVSPARKGLRKpLobIIIAIIIB, getTableauVSPARKGLRKpMidpoint, getTableauVSPARKGLRKpSymmetric, + getTableauVSPARKGLRKpSymmetricLob, getTableauVSPARKGLRKpSymplectic, + getTableauVSPARKLobIIIAIIIB2pLobIIIAIIIB, + getTableauVSPARKLobIIIAIIIB3pLobIIIAIIIB, + getTableauVSPARKLobIIIAIIIB4pLobIIIAIIIB, getTableauVSPARKLobIIIAIIIB2pSymmetric, getTableauVSPARKLobIIIAIIIB3pSymmetric, - getTableauVSPARKLobIIIAIIIB4pSymmetric + getTableauVSPARKLobIIIAIIIB4pSymmetric, + getTableauVSPARKLobIIIAIIIB2pSymmetricLob, + getTableauVSPARKLobIIIAIIIB3pSymmetricLob, + getTableauVSPARKLobIIIAIIIB4pSymmetricLob include("tableaus/tableaus_vspark_primary.jl") diff --git a/src/tableaus/tableaus_vspark_primary.jl b/src/tableaus/tableaus_vspark_primary.jl index ed274e208..7b5d62b54 100644 --- a/src/tableaus/tableaus_vspark_primary.jl +++ b/src/tableaus/tableaus_vspark_primary.jl @@ -126,6 +126,189 @@ function getTableauVSPARKGLRKpSymmetric(s) end + +function getTableauVSPARKSymmetricLobProjection(name, q::CoefficientsRK{T}, p::CoefficientsRK{T}, d=[]; R∞=1) where {T} + + @assert q.s == p.s + + o = min(q.o, p.o) + + loba = getCoefficientsLobIIIA2() + lobb = getCoefficientsLobIIIB2() + + a_q = q.a + a_p = p.a + b_q = q.b + b_p = p.b + + a_q̃ = vcat(zero(q.b)', q.b') + a_p̃ = vcat(zero(p.b)', p.b') +# a_p̃ = vcat(p.b' ./ 2, p.b' ./ 2) + + β_q = [0.5, R∞*0.5] + β_p = [0.5, R∞*0.5] + + α_q̃ = loba.a + α_p̃ = lobb.a + + α_q = zeros(T, q.s, 2) + α_p = zeros(T, p.s, 2) + + for i in 1:q.s + for j in 1:2 + α_q[i,j] = β_q[j] / b_p[i] * ( b_p[i] - a_p̃[j,i] ) + α_p[i,j] = β_p[j] / b_q[i] * ( b_q[i] - a_q̃[j,i] ) + end + end + + c_q = q.c + c_p = p.c + c_λ = [ 0.0, 1.0] + d_λ = [ 0.5, 0.5] + + ω_λ = T[0 0 1] + δ_λ = T[-1 +1] + + + if length(d) == 0 + return TableauVSPARKprimary(name, o, + a_q, a_p, α_q, α_p, + a_q̃, a_p̃, α_q̃, α_p̃, + b_q, b_p, β_q, β_p, + c_q, c_p, c_λ, d_λ, + ω_λ, δ_λ) + else + @assert length(d) == q.s == p.s + + return TableauVSPARKprimary(name, o, + a_q, a_p, α_q, α_p, + a_q̃, a_p̃, α_q̃, α_p̃, + b_q, b_p, β_q, β_p, + c_q, c_p, c_λ, d_λ, + ω_λ, δ_λ, d) + end + +end + + +"Tableau for Gauss-Lobatto IIIA-IIIB method with two stages and symmetric projection." +function getTableauVSPARKLobIIIAIIIB2pSymmetricLob() + getTableauVSPARKSymmetricLobProjection(:LobIIIAIIIB2pSymmetricLob, getCoefficientsLobIIIA2(), getCoefficientsLobIIIB2(), get_lobatto_d_vector(2); R∞=-1) +end + +"Tableau for Gauss-Lobatto IIIA-IIIB method with three stages and symmetric projection." +function getTableauVSPARKLobIIIAIIIB3pSymmetricLob() + getTableauVSPARKSymmetricLobProjection(:LobIIIAIIIB3pSymmetricLob, getCoefficientsLobIIIA3(), getCoefficientsLobIIIB3(), get_lobatto_d_vector(3); R∞=+1) +end + +"Tableau for Gauss-Lobatto IIIA-IIIB method with four stages and symmetric projection." +function getTableauVSPARKLobIIIAIIIB4pSymmetricLob() + getTableauVSPARKSymmetricLobProjection(:LobIIIAIIIB4pSymmetricLob, getCoefficientsLobIIIA4(), getCoefficientsLobIIIB4(), get_lobatto_d_vector(4); R∞=-1) +end + +"Tableau for Gauss-Legendre method with s stages and symplectic projection." +function getTableauVSPARKGLRKpSymmetricLob(s) + glrk = getCoefficientsGLRK(s) + getTableauVSPARKSymmetricLobProjection(Symbol("vpglrk", s, "pSymmetricLob"), glrk, glrk; R∞=(-1)^s) +end + + + +function getTableauVSPARKLobIIIAIIIBProjection(name, q::CoefficientsRK{T}, p::CoefficientsRK{T}, d=[]; R∞=1) where {T} + + @assert q.s == p.s + + s = q.s + o = min(q.o, p.o) + + loba = getCoefficientsLobIIIA2() + lobb = getCoefficientsLobIIIB2() + + a_q = q.a + a_p = p.a + b_q = q.b + b_p = p.b + + # α_q = hcat(ones(s)/2, zeros(s)) + # α_p = hcat(ones(s)/2, zeros(s)) + + α_q = zeros(T, s, 2) + α_p = zeros(T, s, 2) + for i in 1:s + α_q[i,:] .= loba.b + α_p[i,:] .= lobb.b + end + + β_q = loba.b + β_p = lobb.b + + α_q̃ = loba.a + α_p̃ = lobb.a + + a_q̃ = zeros(T, 2, s) + a_p̃ = zeros(T, 2, s) + + for i in 1:2 + for j in 1:s + a_q̃[i,j] = b_q[j] / β_p[i] * ( β_p[i] - α_p[j,i] ) + a_p̃[i,j] = b_p[j] / β_q[i] * ( β_q[i] - α_q[j,i] ) + end + end + + c_q = q.c + c_p = p.c + c_λ = loba.c + d_λ = loba.b + + ω_λ = T[0.5 0.5 0.0 + 0.0 0.0 1.0] + δ_λ = zeros(T, 0, 1) + + + if length(d) == 0 + return TableauVSPARKprimary(name, o, + a_q, a_p, α_q, α_p, + a_q̃, a_p̃, α_q̃, α_p̃, + b_q, b_p, β_q, β_p, + c_q, c_p, c_λ, d_λ, + ω_λ, δ_λ) + else + @assert length(d) == q.s == p.s + + return TableauVSPARKprimary(name, o, + a_q, a_p, α_q, α_p, + a_q̃, a_p̃, α_q̃, α_p̃, + b_q, b_p, β_q, β_p, + c_q, c_p, c_λ, d_λ, + ω_λ, δ_λ, d) + end + +end + + +"Tableau for Gauss-Lobatto IIIA-IIIB method with two stages and symmetric projection." +function getTableauVSPARKLobIIIAIIIB2pLobIIIAIIIB() + getTableauVSPARKSymmetricLobProjection(:LobIIIAIIIB2pSymmetricLob, getCoefficientsLobIIIA2(), getCoefficientsLobIIIB2(), get_lobatto_d_vector(2); R∞=-1) +end + +"Tableau for Gauss-Lobatto IIIA-IIIB method with three stages and symmetric projection." +function getTableauVSPARKLobIIIAIIIB3pLobIIIAIIIB() + getTableauVSPARKSymmetricLobProjection(:LobIIIAIIIB3pSymmetricLob, getCoefficientsLobIIIA3(), getCoefficientsLobIIIB3(), get_lobatto_d_vector(3); R∞=+1) +end + +"Tableau for Gauss-Lobatto IIIA-IIIB method with four stages and symmetric projection." +function getTableauVSPARKLobIIIAIIIB4pLobIIIAIIIB() + getTableauVSPARKSymmetricLobProjection(:LobIIIAIIIB4pSymmetricLob, getCoefficientsLobIIIA4(), getCoefficientsLobIIIB4(), get_lobatto_d_vector(4); R∞=-1) +end + +"Tableau for Gauss-Legendre method with s stages and symplectic projection." +function getTableauVSPARKGLRKpLobIIIAIIIB(s) + glrk = getCoefficientsGLRK(s) + getTableauVSPARKLobIIIAIIIBProjection(Symbol("vpglrk", s, "pLobIIIAIIIB"), glrk, glrk; R∞=(-1)^s) +end + + + function getTableauVSPARKSymplecticProjection(name, q::CoefficientsRK{T}, p::CoefficientsRK{T}, d=[]; R∞=+1) where {T} la = getCoefficientsLobIIIA2() From 7ed43b61a521b93e70cc4a292f5116dc9bf6a809 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 14:43:00 +0200 Subject: [PATCH 29/30] Add Gauss-Legendre tableaus for implicit partitioned Runge-Kutta methods. --- docs/src/releasenotes.md | 1 + src/Tableaus.jl | 4 ++++ src/tableaus/tableaus_iprk.jl | 6 ++++++ 3 files changed, 11 insertions(+) create mode 100644 src/tableaus/tableaus_iprk.jl diff --git a/docs/src/releasenotes.md b/docs/src/releasenotes.md index 0ed9be907..fa7bee2ce 100644 --- a/docs/src/releasenotes.md +++ b/docs/src/releasenotes.md @@ -7,6 +7,7 @@ * Atomic solutions can now store a NamedTuple of internal variables of the integrator, including nonlinear solver output * Output of internal variables has been added to VPRK integrators +* Add Gauss-Legendre tableaus for implicit partitioned Runge-Kutta methods ### Fixes diff --git a/src/Tableaus.jl b/src/Tableaus.jl index 83feabab1..fc26ba36a 100644 --- a/src/Tableaus.jl +++ b/src/Tableaus.jl @@ -101,6 +101,10 @@ module Tableaus export getTableauSPARKGLRK, getTableauSPARKLobIIIAIIIB, getTableauSPARKGLRKLobIIIAIIIB + include("tableaus/tableaus_iprk.jl") + + export getTableauIPGLRK + include("tableaus/tableaus_spark.jl") export getTableauVPGLRK, diff --git a/src/tableaus/tableaus_iprk.jl b/src/tableaus/tableaus_iprk.jl new file mode 100644 index 000000000..8bd91bf33 --- /dev/null +++ b/src/tableaus/tableaus_iprk.jl @@ -0,0 +1,6 @@ + +"Gauss-Legendre Runge-Kutta" +function getTableauIPGLRK(s::Int) + TableauIPRK(Symbol("IPGLRK", s), 2^s, TableauFIRK(getCoefficientsGLRK(s))) + +end From dfe9ad36c2ec878eb69146056e9ba6f264332460 Mon Sep 17 00:00:00 2001 From: Michael Kraus Date: Tue, 22 Sep 2020 15:55:58 +0200 Subject: [PATCH 30/30] Fix in IPRK tableaus. --- src/tableaus/tableaus_eprk.jl | 1 - src/tableaus/tableaus_iprk.jl | 3 +-- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/tableaus/tableaus_eprk.jl b/src/tableaus/tableaus_eprk.jl index c706e0b50..8c778cb55 100644 --- a/src/tableaus/tableaus_eprk.jl +++ b/src/tableaus/tableaus_eprk.jl @@ -1,5 +1,4 @@ - function getCoefficientsSymplecticEulerForward() a = [[0.0 0.0] [1.0 0.0]] diff --git a/src/tableaus/tableaus_iprk.jl b/src/tableaus/tableaus_iprk.jl index 8bd91bf33..2c64e33c9 100644 --- a/src/tableaus/tableaus_iprk.jl +++ b/src/tableaus/tableaus_iprk.jl @@ -1,6 +1,5 @@ "Gauss-Legendre Runge-Kutta" function getTableauIPGLRK(s::Int) - TableauIPRK(Symbol("IPGLRK", s), 2^s, TableauFIRK(getCoefficientsGLRK(s))) - + TableauIPRK(Symbol("IPGLRK", s), 2^s, getCoefficientsGLRK(s)) end