diff --git a/src/integration.jl b/src/integration.jl index 2c9d936f..89520c61 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -25,6 +25,7 @@ end """ boltzmann(prob::CauchyProblem) -> DifferentialEquations.ODEProblem + boltzmann(prob::SorptivityProblem) -> DifferentialEquations.ODEProblem Transform `prob` into an ODE problem in terms of the Boltzmann variable `o`. @@ -32,40 +33,35 @@ The ODE problem is set up to terminate automatically (`ReturnCode.Terminated`) w See also: [`DifferentialEquations`](https://diffeq.sciml.ai/stable/) """ -function boltzmann(prob::CauchyProblem) - u0 = @SVector [prob.b, prob.d_dob] +function boltzmann(prob::Union{CauchyProblem, SorptivityProblem}) + if prob isa CauchyProblem + u0 = @SVector [prob.b, prob.d_dob] + elseif prob isa SorptivityProblem + u0 = @SVector [prob.b, d_do(prob.eq, prob.b, prob.S)] + end + ob = float(prob.ob) + settled = DiscreteCallback(let direction = monotonicity(prob) (u, t, integrator) -> direction * u[2] ≤ zero(u[2]) end, terminate!, save_positions = (false, false)) - ODEProblem(boltzmann(prob.eq), u0, (ob, typemax(ob)), callback = settled) -end -function boltzmann(prob::SorptivityProblem) - boltzmann(CauchyProblem(prob.eq, - b = prob.b, - d_dob = d_do(prob.eq, prob.b, prob.S), - ob = prob.ob)) + ODEProblem(boltzmann(prob.eq), u0, (ob, typemax(ob)), callback = settled) end -monotonicity(odeprob::ODEProblem)::Int = sign(odeprob.u0[2]) - const _ODE_ALG = RadauIIA5() const _ODE_MAXITERS = 1000 -function _init(odeprob::ODEProblem, +function _init(prob::Union{CauchyProblem, SorptivityProblem}, ::BoltzmannODE; - i = nothing, - abstol = 0, + limit = nothing, verbose = true) - if !isnothing(i) - direction = monotonicity(odeprob) - - past_limit = DiscreteCallback(let direction = direction, - limit = i + direction * abstol + odeprob = boltzmann(prob) + if !isnothing(limit) + past_limit = DiscreteCallback(let direction = monotonicity(prob), limit = limit (u, t, integrator) -> direction * u[1] > direction * limit end, terminate!, @@ -80,14 +76,6 @@ function _init(odeprob::ODEProblem, return init(odeprob, _ODE_ALG, maxiters = _ODE_MAXITERS, verbose = verbose) end -function _init(prob::Union{CauchyProblem, SorptivityProblem}, - alg::BoltzmannODE; - i = nothing, - abstol = 0, - verbose = true) - _init(boltzmann(prob), alg, i = i, abstol = abstol, verbose = verbose) -end - function _reinit!(integrator, prob::CauchyProblem) @assert sign(prob.d_dob) == sign(integrator.sol.u[1][2]) reinit!(integrator, @SVector [prob.b, prob.d_dob]) diff --git a/src/shooting.jl b/src/shooting.jl index 1a84f3c1..1ac8cb62 100644 --- a/src/shooting.jl +++ b/src/shooting.jl @@ -1,35 +1,3 @@ -function _shoot!(integrator, prob::Union{CauchyProblem, SorptivityProblem}; i, abstol) - direction = monotonicity(prob) - limit = i + direction * abstol - - integrator = _reinit!(integrator, prob) - - solve!(integrator) - - @assert integrator.sol.retcode != ReturnCode.Success - - resid = direction * typemax(i) - - if integrator.sol.retcode == ReturnCode.Terminated && - direction * integrator.sol.u[end][1] <= direction * limit - resid = integrator.sol.u[end][1] - i - end - - return integrator, resid -end - -function _init(prob::DirichletProblem, alg::BoltzmannODE; d_dob, abstol) - return _init(CauchyProblem(prob.eq, b = prob.b, d_dob = d_dob, ob = prob.ob), - alg, - i = prob.i, abstol = abstol, verbose = false) -end - -function _shoot!(integrator, prob::DirichletProblem; d_dob, abstol) - return _shoot!(integrator, - CauchyProblem(prob.eq, b = prob.b, d_dob = d_dob, ob = prob.ob), i = prob.i, - abstol = abstol) -end - """ solve(prob::DirichletProblem[, alg::BoltzmannODE; abstol, maxiters, d_dob_hint]) -> Solution @@ -64,28 +32,45 @@ function solve(prob::DirichletProblem, alg::BoltzmannODE = BoltzmannODE(); d_dob_hint = d_do(prob, :b_hint) end + direction = monotonicity(prob) + limit = prob.i + direction * abstol resid = prob.b - prob.i - integrator = _init(prob, alg, d_dob = d_dob_hint, abstol = abstol) + integrator = _init(CauchyProblem(prob.eq, b = prob.b, d_dob = d_dob_hint, ob = prob.ob), + alg, + limit = limit, + verbose = false) if abs(resid) ≤ abstol + _reinit!(integrator, + CauchyProblem(prob.eq, b = prob.b, d_dob = zero(d_dob_hint), ob = prob.ob)) solve!(integrator) + @assert integrator.sol.retcode != ReturnCode.Success retcode = integrator.sol.retcode == ReturnCode.Terminated ? ReturnCode.Success : integrator.sol.retcode if verbose && !SciMLBase.successful_retcode(integrator.sol) @warn "Problem has a trivial solution but failed to obtain it" end + return Solution(integrator.sol, prob, alg, _retcode = retcode, _niter = 0) end d_dob_trial = bracket_bisect(zero(d_dob_hint), d_dob_hint, resid) for niter in 1:maxiters - integrator, resid = _shoot!(integrator, - prob, - d_dob = d_dob_trial(resid), - abstol = abstol) + _reinit!(integrator, + CauchyProblem(prob.eq, b = prob.b, d_dob = d_dob_trial(resid), ob = prob.ob)) + solve!(integrator) + + @assert integrator.sol.retcode != ReturnCode.Success + if integrator.sol.retcode == ReturnCode.Terminated && + direction * integrator.sol.u[end][1] ≤ direction * limit + resid = integrator.sol.u[end][1] - prob.i + else + resid = direction * typemax(prob.i) + end + if abs(resid) ≤ abstol return Solution(integrator.sol, prob, @@ -105,24 +90,6 @@ function solve(prob::DirichletProblem, alg::BoltzmannODE = BoltzmannODE(); _niter = maxiters) end -function _init(prob::FlowrateProblem, alg::BoltzmannODE; b, abstol) - ob = !iszero(prob.ob) ? prob.ob : 1e-6 - - return _init(SorptivityProblem(prob.eq, b = b, S = 2prob.Qb / prob._αh / ob, ob = ob), - alg, - i = prob.i, - abstol = abstol, - verbose = false) -end - -function _shoot!(integrator, prob::FlowrateProblem; b, abstol) - ob = !iszero(prob.ob) ? prob.ob : 1e-6 - - return _shoot!(integrator, - SorptivityProblem(prob.eq, b = b, S = 2prob.Qb / prob._αh / ob, ob = ob), - i = prob.i, abstol = abstol) -end - """ solve(prob::FlowrateProblem[, BoltzmannODE; abstol, maxiters, b_hint]) -> Solution @@ -157,19 +124,29 @@ function solve(prob::FlowrateProblem; alg::BoltzmannODE = BoltzmannODE(), b_hint = prob.i - oneunit(prob.i) * monotonicity(prob) end + ob = iszero(prob.ob) ? 1e-6 : prob.ob + + direction = monotonicity(prob) + limit = prob.i + direction * abstol resid = prob.i - oneunit(prob.i) * monotonicity(prob) + S = 2prob.Qb / prob._αh / ob + + integrator = _init(SorptivityProblem(prob.eq, b = b_hint, S = S, ob = ob), + alg, + limit = limit, + verbose = false) - integrator = _init(prob, alg, b = b_hint, abstol = abstol) + if iszero(direction) + _reinit!(integrator, SorptivityProblem(prob.eq, b = prob.i, S = zero(S), ob = ob)) + solve!(integrator) - if monotonicity(prob) == 0 - integrator, resid = _shoot!(integrator, prob, b = prob.i, abstol = abstol) - @assert iszero(resid) @assert integrator.sol.retcode != ReturnCode.Success retcode = integrator.sol.retcode == ReturnCode.Terminated ? ReturnCode.Success : integrator.sol.retcode if verbose && !SciMLBase.successful_retcode(integrator.sol) @warn "Problem has a trivial solution but failed to obtain it" end + return Solution(integrator.sol, prob, alg, @@ -180,7 +157,18 @@ function solve(prob::FlowrateProblem; alg::BoltzmannODE = BoltzmannODE(), b_trial = bracket_bisect(prob.i, b_hint) for niter in 1:maxiters - integrator, resid = _shoot!(integrator, prob, b = b_trial(resid), abstol = abstol) + _reinit!(integrator, SorptivityProblem(prob.eq, b = b_trial(resid), S = S, ob = ob)) + solve!(integrator) + + @assert integrator.sol.retcode != ReturnCode.Success + if integrator.sol.retcode == ReturnCode.Terminated && + direction * integrator.sol.u[end][1] ≤ direction * limit + resid = integrator.sol.u[end][1] - prob.i + elseif integrator.sol.retcode != ReturnCode.Terminated && integrator.t == ob + resid = -direction * typemax(prob.i) + else + resid = direction * typemax(prob.i) + end if abs(resid) ≤ abstol return Solution(integrator.sol,