Skip to content

Commit

Permalink
Merge pull request #47 from LilithHafner/lh/format
Browse files Browse the repository at this point in the history 
Run JuliaFormatter.format()
  • Loading branch information
@ChrisRackauckas
ChrisRackauckas authored Feb 15, 2024
2 parents ea761ca + 1e3002d commit 83846a6
Show file tree
Hide file tree
Showing 4 changed files with 101 additions and 99 deletions.
122 changes: 61 additions & 61 deletions src/DASSL.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,22 +23,22 @@ mutable struct JacData
end

function dasslStep(channel,
F,
y0::AbstractVector{T},
tstart::Real;
reltol = 1e-3,
abstol = 1e-5,
initstep = 1e-4,
maxstep = Inf,
minstep = 0,
maxorder = MAXORDER,
dy0 = zero(y0),
tstop = Inf,
norm = dassl_norm,
weights = dassl_weights,
factorize_jacobian = true, # whether to store factorized version of jacobian
jacobian = numerical_jacobian(F, reltol, abstol, weights), # computes the quantity F/dy+a*dF/dy'
args...) where {T <: Number}
F,
y0::AbstractVector{T},
tstart::Real;
reltol = 1e-3,
abstol = 1e-5,
initstep = 1e-4,
maxstep = Inf,
minstep = 0,
maxorder = MAXORDER,
dy0 = zero(y0),
tstop = Inf,
norm = dassl_norm,
weights = dassl_weights,
factorize_jacobian = true, # whether to store factorized version of jacobian
jacobian = numerical_jacobian(F, reltol, abstol, weights), # computes the quantity F/dy+a*dF/dy'
args...) where {T <: Number}
n = length(y0)
jd = JacData(zero(tstart), zeros(T, n, n)) # generate a dummy
# jacobian, it will be replaced by
Expand All @@ -65,8 +65,8 @@ function dasslStep(channel,
# y0. Basically we run a one iteration of stepper and use the
# result as the initial derivative
(_, _, _, dyout[1], _) = stepper(1, tout, yout, dyout, 10 * eps(one(tstart)), F, jd,
computejac, weights(y0, 1, 1),
v -> norm(v, weights(y0, 1, 1)), 1)
computejac, weights(y0, 1, 1),
v -> norm(v, weights(y0, 1, 1)), 1)

while tout[end] < tstop
hmin = max(4 * eps(one(tstart)), minstep)
Expand All @@ -85,7 +85,7 @@ function dasslStep(channel,
normy(v) = norm(v, wt)

(status, err, yn, dyn, jd) = stepper(ord, tout, yout, dyout, h, F, jd, computejac,
wt, normy, maxorder)
wt, normy, maxorder)

if status < 0
# Early failure: Newton iteration failed to converge, reduce
Expand Down Expand Up @@ -180,12 +180,12 @@ function dasslSolve(F, y0::Number, tspan; args...)
end

function newStepOrder(t::AbstractVector,
y::AbstractVector,
normy,
err,
k::Integer,
num_fail::Integer,
maxorder::Integer)
y::AbstractVector,
normy,
err,
k::Integer,
num_fail::Integer,
maxorder::Integer)
if length(t) != length(y)
error("incompatible size of y and t")
end
Expand Down Expand Up @@ -241,11 +241,11 @@ function newStepOrder(t::AbstractVector,
end

function newStepOrderContinuous(t::AbstractVector,
y::AbstractVector,
normy,
err,
k::Integer,
maxorder::Integer)
y::AbstractVector,
normy,
err,
k::Integer,
maxorder::Integer)

# compute the error estimates of methods of order k-2, k-1, k and
# (if possible) k+1
Expand Down Expand Up @@ -334,9 +334,9 @@ end
# here t is an array of times t=[t_1, ..., t_n, t_{n+1}]
# and y is an array of solutions y=[y_1, ..., y_n, y_{n+1}]
function errorEstimates(t::AbstractVector,
y::AbstractVector,
normy,
k::Integer)
y::AbstractVector,
normy,
k::Integer)
nsteps = length(t) # available steps (including counting
# the new n+1'st step)
h = t[end] - t[end - 1] # current step size
Expand All @@ -361,7 +361,7 @@ function errorEstimates(t::AbstractVector,
hn = diff(t[(end - (k + 2)):end])
if all(hn .== hn[1])
maxd = interpolateHighestDerivative(t[(end - (k + 2)):end],
y[(end - (k + 2)):end])
y[(end - (k + 2)):end])
push!(errors, h^(k + 2) * normy(maxd))
end
end
Expand All @@ -378,16 +378,16 @@ end
# jd is a bunch of auxilary data saved between steps (jacobian and last coefficient 'a')
# wt is a vector of weights of the norm
function stepper(ord::Integer,
t::AbstractVector,
y::AbstractVector,
dy::AbstractVector,
h_next::Real,
F,
jd::JacData,
computejac,
wt,
normy,
maxorder::Integer)
t::AbstractVector,
y::AbstractVector,
dy::AbstractVector,
h_next::Real,
F,
jd::JacData,
computejac,
wt,
normy,
maxorder::Integer)
l = length(y[1]) # the number of dependent variables

# @todo this should be the view of the tail of the arrays t and y
Expand Down Expand Up @@ -433,11 +433,11 @@ function stepper(ord::Integer,

# we compute the corrected value "yc", updating the gradient if necessary
(status, yc, jd) = corrector(jd, # old coefficient a and jacobian
a, # current coefficient a
jac_new, # this function is called when new jacobian is needed
y0, # starting point for modified newton
f_newton, # we want to find zeroes of this function
normy) # the norm used to estimate error needs weights
a, # current coefficient a
jac_new, # this function is called when new jacobian is needed
y0, # starting point for modified newton
f_newton, # we want to find zeroes of this function
normy) # the norm used to estimate error needs weights

alpha = Array{eltype(t)}(undef, ord + 1)

Expand Down Expand Up @@ -471,11 +471,11 @@ end
# the jacobian g_old and a_old.

function corrector(jd::JacData,
a_new::T,
jac_new,
y0::AbstractVector{Ty},
f_newton,
normy) where {T, Ty}
a_new::T,
jac_new,
y0::AbstractVector{Ty},
f_newton,
normy) where {T, Ty}

# if jd.a == 0 the new jacobian is always computed, independently
# of the value of a_new
Expand Down Expand Up @@ -510,8 +510,8 @@ end
# set back to y0. Status tells if the fixed point was obtained
# (status==0) or not (status==-1).
function newton_iteration(f,
y0::AbstractVector{T},
normy) where {T <: Number}
y0::AbstractVector{T},
normy) where {T <: Number}

# first guess comes from the predictor method, then we compute the
# second guess to get the norm1
Expand Down Expand Up @@ -571,8 +571,8 @@ end

# returns the value of the interpolation polynomial at the point x0
function interpolateAt(x::AbstractVector{T},
y::AbstractVector,
x0::T) where {T <: Real}
y::AbstractVector,
x0::T) where {T <: Real}
if length(x) != length(y)
error("x and y have to be of the same size.")
end
Expand All @@ -597,8 +597,8 @@ end
# returns the value of the derivative of the interpolation polynomial
# at the point x0
function interpolateDerivativeAt(x::AbstractVector{T},
y::AbstractVector,
x0::T) where {T <: Real}
y::AbstractVector,
x0::T) where {T <: Real}
if length(x) != length(y)
error("x and y have to be of the same size.")
end
Expand Down Expand Up @@ -632,7 +632,7 @@ end
# p(x)=a_{k-1}*x^{k-1}+...a_1*x+a_0 then this function returns the
# k-th derivative of p, i.e. (k-1)!*a_{k-1}
function interpolateHighestDerivative(x::AbstractVector,
y::AbstractVector)
y::AbstractVector)
if length(x) != length(y)
error("x and y have to be of the same size.")
end
Expand Down
24 changes: 12 additions & 12 deletions src/common.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@ end
dassl(; maxorder = 6, factorize_jacobian = true) = dassl(maxorder, factorize_jacobian)

function solve(prob::DiffEqBase.AbstractDAEProblem{uType, duType, tupType, isinplace},
alg::DASSLDAEAlgorithm, args...; timeseries_errors = true,
abstol = 1e-5, reltol = 1e-3, dt = 1e-4, dtmin = 0.0, dtmax = Inf,
callback = nothing, kwargs...) where {uType, duType, tupType, isinplace}
alg::DASSLDAEAlgorithm, args...; timeseries_errors = true,
abstol = 1e-5, reltol = 1e-3, dt = 1e-4, dtmin = 0.0, dtmax = Inf,
callback = nothing, kwargs...) where {uType, duType, tupType, isinplace}
tType = eltype(tupType)

if callback != nothing || :callback in keys(prob.kwargs)
Expand All @@ -32,14 +32,14 @@ function solve(prob::DiffEqBase.AbstractDAEProblem{uType, duType, tupType, isinp
### Finishing Routine

ts, timeseries, dus = dasslSolve(f, prob.u0, tspan,
abstol = abstol,
reltol = reltol,
maxstep = dtmax,
minstep = dtmin,
initstep = dt,
dy0 = prob.du0,
maxorder = alg.maxorder,
factorize_jacobian = alg.factorize_jacobian)
abstol = abstol,
reltol = reltol,
maxstep = dtmax,
minstep = dtmin,
initstep = dt,
dy0 = prob.du0,
maxorder = alg.maxorder,
factorize_jacobian = alg.factorize_jacobian)
#=
timeseries = Vector{uType}(0)
if typeof(prob.u0)<:Number
Expand All @@ -53,5 +53,5 @@ function solve(prob::DiffEqBase.AbstractDAEProblem{uType, duType, tupType, isinp
end
=#
DiffEqBase.build_solution(prob, alg, ts, timeseries, du = dus,
timeseries_errors = timeseries_errors)
timeseries_errors = timeseries_errors)
end
48 changes: 24 additions & 24 deletions test/convergence.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,11 @@ using DASSL
# analytic and numerical solutions. The third element in the returned
# touple is the time it took to obtain the numerical solution.
function dasslTestConvergence(F::Function, # equation to solve
y0::Vector{T}, # initial data
tspan::Vector{T}, # time span of a solution
sol::Function, # analytic solution, for comparison with numerical solution
rtol_range::Vector{T}, # vector of relative tolerances
atol_range::Vector{T}) where {T <: Number} # vector of absolute tolerances
y0::Vector{T}, # initial data
tspan::Vector{T}, # time span of a solution
sol::Function, # analytic solution, for comparison with numerical solution
rtol_range::Vector{T}, # vector of relative tolerances
atol_range::Vector{T}) where {T <: Number} # vector of absolute tolerances
if length(rtol_range) != length(atol_range)
error("The table of relative errors and absolute errors should be of the same size.")
end
Expand All @@ -23,8 +23,8 @@ function dasslTestConvergence(F::Function, # equation to solve

for i in 1:n
times[i] = @elapsed (tn, yn) = dasslSolve(F, y0, tspan,
rtol = rtol_range[i],
atol = atol_range[i])
rtol = rtol_range[i],
atol = atol_range[i])
k = length(tn)
delta_rel = zeros(T, k)
delta_abs = zeros(T, k)
Expand All @@ -43,11 +43,11 @@ function dasslTestConvergence(F::Function, # equation to solve
end

function dasslDrawConvergence(F::Function, # equation to solve
y0::Vector{T}, # initial data
tspan::Vector{T}, # time span of a solution
sol::Function, # analytic solution, for comparison with numerical solution
rtol::Vector{T}, # vector of relative tolerances
atol::Vector{T}) where {T <: Number} # vector of absolute tolerances
y0::Vector{T}, # initial data
tspan::Vector{T}, # time span of a solution
sol::Function, # analytic solution, for comparison with numerical solution
rtol::Vector{T}, # vector of relative tolerances
atol::Vector{T}) where {T <: Number} # vector of absolute tolerances

# run the dasslSolve to trigger the compilation
dasslSolve(F, y0, tspan)
Expand All @@ -58,20 +58,20 @@ function dasslDrawConvergence(F::Function, # equation to solve
tbl = Table(1, 3)

tbl[1, 1] = FramedPlot(title = "Absolute error",
xlabel = "Absolute tolerance",
ylabel = "Absolute error",
xlog = true,
ylog = true)
xlabel = "Absolute tolerance",
ylabel = "Absolute error",
xlog = true,
ylog = true)
tbl[1, 2] = FramedPlot(title = "Relative error",
xlabel = "Relative tolerance",
ylabel = "Relative error",
xlog = true,
ylog = true)
xlabel = "Relative tolerance",
ylabel = "Relative error",
xlog = true,
ylog = true)
tbl[1, 3] = FramedPlot(title = "Execution time",
xlabel = "Relative error",
ylabel = "Elapsed time",
xlog = true,
ylog = true)
xlabel = "Relative error",
ylabel = "Elapsed time",
xlog = true,
ylog = true)

add(tbl[1, 1], Points(atol, aerr))
add(tbl[1, 2], Points(rtol, rerr))
Expand Down
6 changes: 4 additions & 2 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ using DASSL, Test

# analytical jacobian version (vector)
(tna, yna, dyna) = dasslSolve(F, [sol(0.0)], tspan, maxorder = order, Fy = Fy,
Fdy = Fdy)
Fdy = Fdy)
aerror = maximum(abs.(map(first, yn) - sol(tn)))
rerror = maximum(abs.(map(first, yn) - sol(tn)) ./ abs.(sol(tn)))
nsteps = length(tn)
Expand All @@ -40,4 +40,6 @@ using DASSL, Test
end
end

@testset "Testing common interface" begin include("common.jl") end
@testset "Testing common interface" begin
include("common.jl")
end

0 comments on commit 83846a6

Please sign in to comment.