Skip to content
This repository has been archived by the owner on Sep 9, 2024. It is now read-only.

Commit

Permalink
finally got it working
Browse files Browse the repository at this point in the history
  • Loading branch information
@mauro3
mauro3 committed Jun 10, 2015
1 parent c1b35eb commit 1e0fff1
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 64 deletions.
28 changes: 15 additions & 13 deletions examples/bruss1d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,34 +169,36 @@ ic[1:2:2N] = u0(x)
ic[2:2:2N] = v0(x)
tspan = T[0.0, 10.0] # integration interval
tspan0 = T[0.0, 0.001] # integration interval
tspan1 = T[0.0, 0.3] # integration interval

#t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac)

## ode23s
# tout, yout = ode23s(fn_ex, ic, tspan0, jacobian=jac_ex)
# @time tout, yout = ode23s(fn_ex, ic, tspan, jacobian=jac_ex)
# @show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)
tout, yout = ode23s(fn_ex, ic, tspan0, jacobian=jac_ex)
@time tout, yout = ode23s(fn_ex, ic, tspan, jacobian=jac_ex)
@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)

# ## DASSL
# t,ydassl = dasslSolve(fn, ic, tspan0, jacobian=jac)
# @time t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac)
# @show norm(abs(ydassl[end][refsolinds]-refsol)./refsol, Inf)
t,ydassl = dasslSolve(fn, ic, tspan0, jacobian=jac)
@time t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac)
@show norm(abs(ydassl[end][refsolinds]-refsol)./refsol, Inf)

## ROSW

# # No jac works but is slower and gives much higher precision than
# # other methods.
# tout, yout = ode_rosw(fn!, ic, tspan0)
# @time tout, yout = ode_rosw(fn!, ic, tspan)
# @show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)
# No jac works but is slower and gives much higher precision than
# other methods.
tout, yout = ode_rosw(fn!, ic, tspan0)
@time tout, yout = ode_rosw(fn!, ic, tspan)
@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)

# Analytic Jacobian
# Analytic Jacobian, same as previous
tout, yout = ode_rosw(fn!, ic, tspan0, jacobian=(jac!, jac!()))
@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=(jac!, jac!()))
@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)

# Numerical colored Jacobian
tout, yout = ode_rosw(fn!, ic, tspan0, jacobian=jac!())
@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=jac!())
@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=jac!(), abstol=1e-3, reltol=1e-5);
@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)

nothing
104 changes: 53 additions & 51 deletions src/rosw.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,12 +213,12 @@ function oderosw_adapt{N,S}(fn!, x0::AbstractVector, tspan, btab::TableauRosW{N,
k = Array(Tx, S) # stage variables
allocate!(k, x0, dof)
ks = zeros(Exf,dof) # work vector for one k
if isa(jacobian, Tuple) # Jacobian function and sparse storage provided
if isa(jacobian, Tuple) # Jacobian function and sparse pattern/storage provided
jac_store = jacobian[2]
jacobian = jacobian[1]
elseif isa(jacobian, SparseMatrixCSC) # Sparse storage provided, make colored Jacobian (TODO)
jac_store = jacobian
jacobian = numerical_jacobian(fn!, reltol, abstol) #, jac_store, jac_store)
jacobian = numerical_jacobian(fn!, reltol, abstol, jac_store)
else # nothing provided, make dense numerical jacobian
jac_store = zeros(Exf,dof,dof) # Jacobian storage
end
Expand Down Expand Up @@ -352,7 +352,6 @@ function rosw_step!{N,S}(xtrial, g!, gprime!, x, dt, dof, btab::TableauRosW_T{N,
# TODO: add option to do more Newton iterations
k[s] = A_ldiv_B!(jac, k[s]) # TODO: see whether this is impacted by https://github.com/JuliaLang/julia/issues/10787
# and https://github.com/JuliaLang/julia/issues/11325
# k[s][:] = jac_store\k[s]
for d=1:dof
k[s][d] = ks[d]-k[s][d]
end
Expand All @@ -363,7 +362,6 @@ function rosw_step!{N,S}(xtrial, g!, gprime!, x, dt, dof, btab::TableauRosW_T{N,
# k[S][:] = ks - jac\k[S]
# TODO: add option to do more Newton iterations
k[S] = A_ldiv_B!(jac, k[S])
# k[S][:] = jac_store\k[S]
for d=1:dof
k[S][d] = ks[d]-k[S][d]
end
Expand Down Expand Up @@ -415,7 +413,7 @@ function hprime!(res, xi, gprime!, u, udot, btab, dt)
# The res will hold part of the LU factorization. However use the
# returned LU-type.
gprime!(res, u, udot, 1./(dt*btab.γii))
if true # issparse(res)
if issparse(res)
# For issues with sparse https://github.com/JuliaLang/julia/issues/7488
# Basically lufact! destroys res and it cannot be used anymore later.
return lufact(res)
Expand All @@ -428,79 +426,83 @@ end
# generate a function that computes approximate jacobian using forward
# finite differences
function numerical_jacobian(fn!, reltol, abstol)
function numjac!(jac, y, dy, a)
function numjac!(jac, y, ydot, a)
# Numerical Jacobian, finite differences, no coloring.
ep = eps(1e6) # this is the machine epsilon # TODO: better value?
ep = eps(1.0) # this is the machine epsilon # TODO: better value?
h = 1/a # h ~ 1/a
wt = reltol*abs(y).+abstol
# Delta for approximation of Jacobian. I removed the
# sign(h_next*dy0) from the definition of delta because it was
# causing trouble when dy0==0 (which happens for ord==1)
edelta = max(abs(y),abs(h*dy),wt)*sqrt(ep)

tmp = similar(y)
tmp1 = similar(y)
tmp2 = similar(y)
# sign(h_next*ydot0) from the definition of delta because it was
# causing trouble when ydot0==0 (which happens for ord==1)
edelta = max(abs(y),abs(h*ydot),wt)*sqrt(ep)

f0 = similar(y)
fn!(f0, y, dy)
fn!(f0, y, ydot)

dy = deepcopy(y)
dydot = deepcopy(ydot)
f1 = similar(y)
for i=1:length(y)
copy!(tmp1, y)
copy!(tmp2, dy)
tmp1[i] += edelta[i]
tmp2[i] += a*edelta[i]
fn!(tmp, tmp1, tmp2)
if false #isa(jac, SparseMatrixCSC)
for nz in nzrange(jac,i)
nonzeros(jac)[nz] = tmpy[rowvals(jac)[nz]]
end
else
for j=1:length(y)
jac[j,i] = (tmp[j]-f0[j])/edelta[i]
end
dy[i] += edelta[i]
dydot[i] += a*edelta[i]
fn!(f1, dy, dydot)
for j=1:length(y)
jac[j,i] = (f1[j]-f0[j])/edelta[i]
end
# remove delta again
dy[i] -= edelta[i]
dydot[i] -= a*edelta[i]
end
return nothing
end
end

function numerical_jacobian(fn!, reltol, abstol, JPattern1::SparseMatrixCSC, JPattern2::SparseMatrixCSC)
# JPattern1: sparsity pattern of dfn/dy
# JPattern2: sparsity pattern of dfn/dydot

# S1 = color_jac(JPattern1) # seed matrix
# S2 = color_jac(JPattern2) # seed matrix
# nc1 = size(S1,2) # number of colors
# nc2 = size(S2,2) # number of colors
if VERSION < v"0.4.0-dev"
# TODO: remove these two for 0.4
rowvals(S::SparseMatrixCSC) = S.rowval
nzrange(S::SparseMatrixCSC, col::Integer) = S.colptr[col]:(S.colptr[col+1]-1)
end
function numerical_jacobian(fn!, reltol, abstol, JPattern::SparseMatrixCSC)
# JPattern: sparsity pattern of dfn/dy + dfn/dydot
# Typically this matrix will be used for storage of the
# Jacobian as well.

# TODO: there are probably more clever ways to do this:
S = color_jac(JPattern1+JPattern2)
S = color_jac(JPattern)
nc = size(S,2)
function numjac_c!(jac, y, ydot, a)
# updates jac in-place
ep = eps(1e6) # this is the machine epsilon # TODO: better value?
ep = eps(1.0) # this is the machine epsilon # TODO: better value?
h = 1/a # h ~ 1/a
wt = reltol*abs(y).+abstol

edelta = max(abs(y),abs(h*ydot),wt)*sqrt(ep)
a_edelta = a*edelta

f0 = similar(y)
fn!(f0, y, ydot)

dof = length(y)
res0 = similar(y)
fn!(res0, y, ydot)
dy = deepcopy(y) # y+dy
dydot = deepcopy(y) # y+dydot
res1 = similar(y) # fn(y+dy, ydot)
res2 = similar(y) # fn(y, ydot+dydot)
dydot = deepcopy(ydot) # y+dydot
f1 = similar(y) # fn(y+dy, ydot)
for c in 1:nc
add_delta!(dy, S, c, edelta) # adds a delta to all indices of color c
add_delta!(dydot, S, c, edelta) # adds a delta to all indices of color c
fn!(res1, dy, ydot)
fn!(res2, y, dydot)
for j=1:dof
res1[j] = ( res1[j]-res0[j] + a*(res2[j]-res0[j]) )/edelta[j]
add_delta!(dydot, S, c, a_edelta) # adds a delta to all indices of color c
fn!(f1, dy, dydot)
for j=1:length(y)
f1[j] = f1[j]-f0[j]
end
recover_jac!(jac, res1, S, c) # recovers the columns of jac which correspond c
recover_jac!(jac, f1, S, c) # recovers the columns of jac which correspond c
remove_delta!(dy, S, c, edelta) # remove delta again
remove_delta!(dydot, S, c, edelta) # remove delta again
remove_delta!(dydot, S, c, a_edelta) # remove delta again
end
# do the division by edelta
nz = nonzeros(jac)
for j =1:size(jac,2)
ed = edelta[j]
for i in nzrange(jac,j)
nz[i] /= ed
end
end
return nothing
end
Expand Down

0 comments on commit 1e0fff1

Please sign in to comment.