Skip to content

Commit

Permalink
Update BoltzmannODE
Browse files Browse the repository at this point in the history
  • Loading branch information
@gerlero
gerlero committed Dec 18, 2023
1 parent d54e22f commit cdd7226
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 87 deletions.
42 changes: 15 additions & 27 deletions src/integration.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,47 +25,43 @@ end

"""
boltzmann(prob::CauchyProblem) -> DifferentialEquations.ODEProblem
boltzmann(prob::SorptivityProblem) -> DifferentialEquations.ODEProblem
Transform `prob` into an ODE problem in terms of the Boltzmann variable `o`.
The ODE problem is set up to terminate automatically (`ReturnCode.Terminated`) when the steady state is reached.
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!,
Expand All @@ -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])
Expand Down
108 changes: 48 additions & 60 deletions src/shooting.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down

0 comments on commit cdd7226

Please sign in to comment.