Add SorptivityProblem

docs/src/problems.md

@@ -8,6 +8,7 @@ CurrentModule = Fronts
 Problem
 DirichletProblem
 FlowrateProblem
+SorptivityProblem
 CauchyProblem
 SorptivityCauchyProblem
 FiniteProblem

src/problems.jl

@@ -158,6 +158,67 @@ end
 
 monotonicity(prob::FlowrateProblem)::Int = -sign(prob.Qb)
 
+"""
+    SorptivityProblem(eq; i, S[, ob]) <: Problem{typeof(eq)}
+    SorptivityProblem(D; i, S[, ob]) <: Problem{typeof(eq)}
+
+Semi-infinite problem with a known initial condition and soprtivity.
+
+# Arguments
+- `eq::Equation`: governing equation.
+- `D`: diffusivity function. Shortcut for `SorptivityProblem(DiffusionEquation(D), ...)`.
+
+# Keyword arguments
+- `i`: initial value.
+- `S`: prescribed sorptivity.
+- `ob=0`: boundary constant for an optional moving boundary. At time `t`, the boundary is located at `ob*√t`. Must be positive if `eq` is radial.
+
+# Examples
+```jldoctest; setup = :(using Fronts)
+julia> D(θ) = θ^4
+D (generic function with 1 method)
+
+julia> prob = Fronts.SorptivityProblem(D, i=0, S=1)
+⎧ ∂θ/∂t = ∂(D(θ)*∂θ/∂r)/∂r, r>0,t>0
+⎨ θ(r,0) = 0, r>0
+⎩ S = 1
+```
+"""
+struct SorptivityProblem{Teq, _T, _To, _TS} <: Problem{Teq}
+    eq::Teq
+    i::_T
+    S::_TS
+    ob::_To
+    function SorptivityProblem(eq::Equation{1}; i, S, ob = 0)
+        new{typeof(eq), typeof(i), typeof(ob), typeof(S)}(eq, i, S, ob)
+    end
+    function SorptivityProblem(eq::Equation; i, S, ob)
+        @argcheck ob > zero(ob)
+        new{typeof(eq), typeof(i), typeof(ob), typeof(S)}(eq, i, S, ob)
+    end
+end
+
+function SorptivityProblem(D; i, S, ob = 0)
+    SorptivityProblem(DiffusionEquation(D), i = i, S = S, ob = ob)
+end
+
+function Base.show(io::IO, prob::SorptivityProblem)
+    if iszero(prob.ob)
+        println(io, "⎧ ", prob.eq, ", r>0,t>0")
+    else
+        println(io, "⎧ ", prob.eq, ", r>rb(t),t>0")
+    end
+    println(io, "⎨ ", prob.eq.sym, "(r,0) = ", prob.i, ", r>0")
+    println(io, "⎩ S = ", prob.S)
+    if !iszero(prob.ob)
+        print(io, "with rb(t) = ", prob.ob, "*√t")
+    end
+end
+
+monotonicity(prob::SorptivityProblem)::Int = -sign(prob.S)
+
+sorptivity(prob::SorptivityProblem) = prob.S
+
 """
     CauchyProblem(eq; b, d_dob[, ob]) <: Problem{typeof(eq)}
     CauchyProblem(D; b, d_dob[, ob]) <: Problem{DiffusionEquation{1}}

src/shooting.jl

@@ -91,7 +91,8 @@ function solve(prob::DirichletProblem, alg::BoltzmannODE = BoltzmannODE();
 end
 
 """
-    solve(prob::FlowrateProblem[, BoltzmannODE; abstol, maxiters, b_hint]) -> Solution
+    solve(prob::FlowrateProblem[, alg::BoltzmannODE; abstol, maxiters, b_hint]) -> Solution
+    solve(prob::SorptivityProblem[, alg::BoltzmannODE; abstol, maxiters, b_hint, verbose]) -> Solution
 
 Solve the problem `prob`.
 
@@ -108,9 +109,10 @@ Solve the problem `prob`.
 GERLERO, G. S.; BERLI, C. L. A.; KLER, P. A. Open-source high-performance software packages for direct and inverse solving of horizontal capillary flow.
 Capillarity, 2023, vol. 6, no. 2, p. 31-40.
 
-See also: [`Solution`](@ref), [`BoltzmannODE`](@ref)
+See also: [`Solution`](@ref), [`BoltzmannODE`](@ref), [`sorptivity`](@ref)
 """
-function solve(prob::FlowrateProblem; alg::BoltzmannODE = BoltzmannODE(),
+function solve(prob::Union{FlowrateProblem, SorptivityProblem},
+        alg::BoltzmannODE = BoltzmannODE();
         abstol = 1e-3,
         maxiters = 100,
         verbose = true)
@@ -124,12 +126,13 @@ 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
+    ob = prob isa FlowrateProblem && 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
+
+    S = prob isa FlowrateProblem ? 2prob.Qb / prob._αh / ob : prob.S
 
     integrator = _init(SorptivityCauchyProblem(prob.eq, b = b_hint, S = S, ob = ob),
         alg,