Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add defect control adaptivity #89

Merged
merged 30 commits into from
Jul 24, 2023
Merged

Add defect control adaptivity #89

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
e2fc944
Add defect control adaptivity
ErikQQY May 29, 2023
33b221e
Fix comments
ErikQQY May 31, 2023
f2dd8f6
format
ErikQQY May 31, 2023
7d028ba
Fix errors
ErikQQY Jun 4, 2023
adcc35a
Fix more errors
ErikQQY Jun 5, 2023
8863e45
Make mesh refinement work
ErikQQY Jun 7, 2023
9f9e437
format
ErikQQY Jun 7, 2023
1bfc50e
Ready for review
ErikQQY Jun 11, 2023
642e5c8
Merge branch 'master' of https://github.com/ErikQQY/BoundaryValueDiff…
ErikQQY Jun 11, 2023
13cbeae
format
ErikQQY Jun 11, 2023
d79c04b
Fix k_interp
ErikQQY Jun 14, 2023
ebfd30e
Merge branch 'master' into master
ErikQQY Jun 15, 2023
6429528
Fix MIRKTableau
ErikQQY Jun 21, 2023
cb7ff37
Done mesh refinement
ErikQQY Jun 22, 2023
3cbfecf
Merge branch 'master' into master
ErikQQY Jun 24, 2023
1f59a70
Refactor
ErikQQY Jun 27, 2023
38d4eaa
Fix sum_stages
ErikQQY Jun 28, 2023
7fece70
Refactor MIRKGeneralCache
ErikQQY Jun 30, 2023
2409831
Refactor k_interp
ErikQQY Jul 1, 2023
89148bb
Fix BVPSystem in loop
ErikQQY Jul 3, 2023
711acd0
x in tableau should be s x s
ErikQQY Jul 4, 2023
2cf2a9f
Fix collocation error
ErikQQY Jul 6, 2023
8d22794
Fix MIRK6 method and skip MIRK tests for now
ErikQQY Jul 6, 2023
7c2cc23
Fianlly working
ErikQQY Jul 7, 2023
dbebe2f
Refactor
ErikQQY Jul 17, 2023
a13f1d4
Refactor half_mesh
ErikQQY Jul 18, 2023
a02c6e0
Refactor
ErikQQY Jul 21, 2023
14d5a34
Fix retcode of NonlinearSolve
ErikQQY Jul 22, 2023
1009899
Merge branch 'master' into master
ErikQQY Jul 23, 2023
6a4860f
All done
ErikQQY Jul 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,5 @@
*.jl.*.cov
*.jl.mem
Manifest.toml
.vscode
.vscode/*
7 changes: 5 additions & 2 deletions src/cache.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,17 @@ function alg_cache{T,U}(alg::MIRK4, S::BVPSystem{T,U})
end
=#

struct MIRK4GeneralCache{kType} <: GeneralMIRKCache
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need this to be mutable? It seems that the array elements are being updated later.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

for i in 1:(N - 1)
    h = x[i + 1] - x[i]
    # Update K
    update_K!(S, cache, TU, i, h)
    for r in 1:s
        x_new = x[i] + c[r] * h
        y_new = (1 - v[r]) * y[i] + v[r] * y[i + 1]
        if r > 1
            y_new += h * sum(j -> X[r, j] * K[j], 1:(r - 1))
        end
        temp = zeros(Float64, M)
        fun!(temp, y_new, S.p, x_new)
        K[r] = temp[:]
    end
    # Update residual
    residual[i] = y[i + 1] - y[i] - h * sum(j -> b[j] * K[j], 1:s)
    cache.k_discrete[i, :] = K[:]
end

The MIRK cache stores the discrete stages from evaluating PHI in the nonlinear solving process. Since we only need discrete stages from the final nonlinear solving step, I think the previous steps can be omitted, but I didn't figure out an executable way, e.g. judging if this step is the last time we evaluate PHI🤔.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's mutating the values but not the cache itself

mutable struct MIRK4GeneralCache{kType} <: GeneralMIRKCache
K::kType
#k_discrete stores discrete stages for each subinterval,
#hence the size of k_discrete is n x (len*s)
k_discrete::AbstractArray
end

@truncate_stacktrace MIRK4GeneralCache

function alg_cache(alg::Union{GeneralMIRK4, MIRK4}, S::BVPSystem{T, U}) where {T, U}
MIRK4GeneralCache([similar(S.y[1]) for i in 1:4])
MIRK4GeneralCache([similar(S.y[1]) for i in 1:4], zeros(S.N, S.order * S.M))
end

struct MIRK6GeneralCache{kType} <: GeneralMIRKCache
Expand Down
92 changes: 9 additions & 83 deletions src/collocation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,96 +36,22 @@ end
S.bc!(S.residual[end], (S.y[1], S.y[end]), S.p, (S.x[1], S.x[end]))
end

#=
@inline function banded_update_K!(S::BVPSystem, cache::AbstractMIRKCache, TU::MIRKTableau, i, h)
M, N, residual, x, y, fun!, order, y_new = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order, S.tmp
c, v, b, X = TU.c, TU.v, TU.b, TU.x
## K, LJ, RJ, Jacobian = cache.K, cache.LJ, cache.RJ, cache.Jacobian
K = cache.K

## index = (M*(i-1)+2):(M*i+1)
## Lindex = (M*(i-1)+1):(M*i)
## Rindex = (M*(i-1)+1+M):(M*i+M)

function Kᵣ!(Kr, y, y₁, r)
x_new = x[i] + c[r]*h
y_new = (1-v[r])*y + v[r]*y₁
if r > 1
y_new += h * sum(j->X[r, j]*K[j], 1:r-1)
end
fun!(x_new, y_new, Kr)
end

# L is the left strip, and R is the right strip
# Lᵢ = -I - hᵢ*Σᵣbᵣ*(∂Kᵣ/∂yᵢ)
# Rᵢ = I - hᵢ*Σᵣbᵣ*(∂Kᵣ/∂y_{i+1})
# From the paper "A Runge-Kutta Type Boundary Value ODE Solver with Defect Control"
# by W.H. Enright and Paul Muir
for r in 1:order
Kᵣ!(K[r], y[i], y[i+1], r)
# ∂Kᵣ/∂yᵢ
### ForwardDiff.jacobian!(LJ[r], (Kr, y₀)->Kᵣ!(Kr, y₀, y[i+1], r), K[r], y[i])
# ∂Kᵣ/∂y_{i+1}
## ForwardDiff.jacobian!(RJ[r], (Kr, y₁)->Kᵣ!(Kr, y[i], y₁, r), K[r], y[i+1])
# h*bᵣ*(∂Kᵣ/∂yᵢ)
## scale!(-b[r]*h, LJ[r])
# hᵢ*Σᵣbᵣ*(∂Kᵣ/∂y_{i+1})
## scale!(-b[r]*h, RJ[r])
# sum them up
## Jacobian[index,Lindex] += LJ[r]
## Jacobian[index,Rindex] += RJ[r]
# fun_jac!(LJ[r], fun!, x_new, y_new, K[r])
# fun_jac!(RJ[r], fun!, x_new, y_new, K[r])
end
# Lᵢ = -I - ...
# Rᵢ = I - ...
## Jacobian[index,Lindex] -= I
## Jacobian[index,Rindex] += I
end

function banded_Φ!(S::BVPSystem, TU::MIRKTableau, cache::AbstractMIRKCache)
order, residual, N, y, x = S.order, S.residual, S.N, S.y, S.x
K, b = cache.K, TU.b
for i in 1:N-1
h = x[i+1] - x[i]
banded_update_K!(S, cache, TU, i, h)
# Update residual
residual[i] = y[i+1] - y[i] - h * sum(j->b[j]*K[j], 1:order)
end
eval_bc_residual!(S)
#ForwardDiff.jacobian!(@view(cache.Jacobian[1:S.M, 1:S.M]), (x,y)->S.bc!(x,y,S.y[1]), residual[1], S.y[1])
#ForwardDiff.jacobian!(@view(cache.Jacobian[(end-S.M+1):end, (end-S.M+1):end]), (x,y)->S.bc!(x,S.y[end],y), residual[end], S.y[end])
#display(cache.Jacobian)
end
=#

@inline function update_K!(S::BVPSystem, cache::AbstractMIRKCache, TU::MIRKTableau, i, h)
M, N, residual, x, y, fun!, order = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order
K, b = cache.K, TU.b
c, v, X = TU.c, TU.v, TU.x

function Kᵣ!(Kr, y, y₁, r)
x_new = x[i] + c[r] * h
y_new = (1 - v[r]) * y + v[r] * y₁
if r > 1
y_new += h * sum(j -> X[r, j] * K[j], 1:(r - 1))
end
fun!(Kr, y_new, S.p, x_new)
end

for r in 1:order
Kᵣ!(K[r], y[i], y[i + 1], r)
end
end

function Φ!(S::BVPSystem{T}, TU::MIRKTableau, cache::AbstractMIRKCache) where {T}
M, N, residual, x, y, fun!, order = S.M, S.N, S.residual, S.x, S.y, S.fun!, S.order
K, b = cache.K, TU.b
c, v, X = TU.c, TU.v, TU.x
for i in 1:(N - 1)
h = x[i + 1] - x[i]
# Update K
update_K!(S, cache, TU, i, h)
for r in 1:order
x_new = x[i] + c[r] * h
y_new = (1 - v[r]) * y[i] + v[r] * y[i + 1]
if r > 1
y_new += h * sum(j -> X[r, j] * K[j], 1:(r - 1))
end
fun!(K[r], y_new, S.p, x_new)
cache.k_discrete[i, ((r - 1) * M + 1):((r - 1) * M + M)] = K[r]
end
# Update residual
residual[i] = y[i + 1] - y[i] - h * sum(j -> b[j] * K[j], 1:order)
end
Expand Down
Loading
Loading