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

Apply W-transform to Rosenbrock23 #2307

Merged
merged 7 commits into from
Sep 12, 2024
Merged

Apply W-transform to Rosenbrock23 #2307

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1438a63
Reapply "Merge remote-tracking branch 'origin/master'"
ChrisRackauckas Aug 4, 2024
adbb347
start fixing addsteps
oscardssmith Aug 7, 2024
1938d52
fix addsteps
oscardssmith Aug 13, 2024
5b03313
fix
oscardssmith Aug 13, 2024
34750c3
Update lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
ChrisRackauckas Aug 16, 2024
0f6931d
Merge branch 'master' into wtransform
oscardssmith Sep 11, 2024
b3d0d4a
build correct jac_config
oscardssmith Sep 12, 2024
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
4 changes: 2 additions & 2 deletions lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -146,7 +146,7 @@ function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits},
assumptions = LinearSolve.OperatorAssumptions(true))

grad_config = build_grad_config(alg, f, tf, du1, t)
jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2, Val(false))
jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
algebraic_vars = f.mass_matrix === I ? nothing :
[all(iszero, x) for x in eachcol(f.mass_matrix)]

Expand Down Expand Up @@ -191,7 +191,7 @@ function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits},
Pl = Pl, Pr = Pr,
assumptions = LinearSolve.OperatorAssumptions(true))
grad_config = build_grad_config(alg, f, tf, du1, t)
jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2, Val(false))
jac_config = build_jac_config(alg, f, uf, du1, uprev, u, tmp, du2)
algebraic_vars = f.mass_matrix === I ? nothing :
[all(iszero, x) for x in eachcol(f.mass_matrix)]

Expand Down
126 changes: 64 additions & 62 deletions lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_perform_step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,8 @@ end
mass_matrix = integrator.f.mass_matrix

# Precalculations
γ = dt * d
dtγ = dt * d
neginvdtγ = -inv(dtγ)
dto2 = dt / 2
dto6 = dt / 6

Expand All @@ -42,7 +43,7 @@ end
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
end

calc_rosenbrock_differentiation!(integrator, cache, γ, γ, repeat_step, false)
calc_rosenbrock_differentiation!(integrator, cache, dtγ, dtγ, repeat_step, true)

calculate_residuals!(weight, fill!(weight, one(eltype(u))), uprev, uprev,
integrator.opts.abstol, integrator.opts.reltol,
Expand All @@ -52,20 +53,20 @@ end
linres = dolinsolve(
integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
solverdata = (; gamma = γ))
solverdata = (; gamma = dtγ))
else
linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
solverdata = (; gamma = γ))
solverdata = (; gamma = dtγ))
end

vecu = _vec(linres.u)
veck₁ = _vec(k₁)

@.. broadcast=false veck₁=-vecu
@.. veck₁ = vecu * neginvdtγ
integrator.stats.nsolve += 1

@.. broadcast=false u=uprev + dto2 * k₁
@.. u=uprev + dto2 * k₁
stage_limiter!(u, integrator, p, t + dto2)
f(f₁, u, p, t + dto2)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
Expand All @@ -76,17 +77,16 @@ end
mul!(_vec(tmp), mass_matrix, _vec(k₁))
end

@.. broadcast=false linsolve_tmp=f₁ - tmp
@.. linsolve_tmp = f₁ - tmp

linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
vecu = _vec(linres.u)
veck2 = _vec(k₂)
veck₂ = _vec(k₂)

@.. broadcast=false veck2=-vecu
@.. veck₂ = vecu * neginvdtγ + veck₁
integrator.stats.nsolve += 1

@.. broadcast=false k₂+=k₁
@.. broadcast=false u=uprev + dt * k₂
@.. u = uprev + dt * k₂
stage_limiter!(u, integrator, p, t + dt)
step_limiter!(u, integrator, p, t + dt)

Expand All @@ -107,7 +107,7 @@ end
linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
vecu = _vec(linres.u)
veck3 = _vec(k₃)
@.. broadcast=false veck3=-vecu
@.. veck3 = vecu * neginvdtγ

integrator.stats.nsolve += 1

Expand All @@ -127,8 +127,8 @@ end

if mass_matrix !== I
algvar = reshape(cache.algebraic_vars, size(u))
@.. broadcast=false atmp=ifelse(algvar, fsallast, false) /
integrator.opts.abstol
invatol = inv(integrator.opts.abstol)
@.. atmp = ifelse(algvar, fsallast, false) * invatol
integrator.EEst += integrator.opts.internalnorm(atmp, t)
end
end
Expand All @@ -145,7 +145,8 @@ end
mass_matrix = integrator.f.mass_matrix

# Precalculations
γ = dt * d
dtγ = dt * d
neginvdtγ = -inv(dtγ)
dto2 = dt / 2
dto6 = dt / 6

Expand All @@ -154,7 +155,7 @@ end
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
end

calc_rosenbrock_differentiation!(integrator, cache, γ, γ, repeat_step, false)
calc_rosenbrock_differentiation!(integrator, cache, dtγ, dtγ, repeat_step, true)

calculate_residuals!(weight, fill!(weight, one(eltype(u))), uprev, uprev,
integrator.opts.abstol, integrator.opts.reltol,
Expand All @@ -164,17 +165,17 @@ end
linres = dolinsolve(
integrator, cache.linsolve; A = nothing, b = _vec(linsolve_tmp),
du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
solverdata = (; gamma = γ))
solverdata = (; gamma = dtγ))
else
linres = dolinsolve(integrator, cache.linsolve; A = W, b = _vec(linsolve_tmp),
du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight,
solverdata = (; gamma = γ))
solverdata = (; gamma = dtγ))
end

vecu = _vec(linres.u)
veck₁ = _vec(k₁)

@.. broadcast=false veck₁=-vecu
@.. veck₁ = vecu * neginvdtγ
integrator.stats.nsolve += 1

@.. broadcast=false u=uprev + dto2 * k₁
Expand All @@ -192,13 +193,12 @@ end

linres = dolinsolve(integrator, linres.cache; b = _vec(linsolve_tmp))
vecu = _vec(linres.u)
veck2 = _vec(k₂)
veck₂ = _vec(k₂)

@.. broadcast=false veck2=-vecu
@.. veck₂ = vecu * neginvdtγ + veck₁
integrator.stats.nsolve += 1

@.. broadcast=false k₂+=k₁
@.. broadcast=false tmp=uprev + dt * k₂
@.. tmp = uprev + dt * k₂
stage_limiter!(u, integrator, p, t + dt)
f(fsallast, tmp, p, t + dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
Expand All @@ -216,7 +216,7 @@ end
vecu = _vec(linres.u)
veck3 = _vec(k₃)

@.. broadcast=false veck3=-vecu
@.. veck3 = vecu * neginvdtγ
integrator.stats.nsolve += 1

@.. broadcast=false u=uprev + dto6 * (k₁ + 4k₂ + k₃)
Expand All @@ -230,8 +230,8 @@ end
integrator.EEst = integrator.opts.internalnorm(atmp, t)

if mass_matrix !== I
@.. broadcast=false atmp=ifelse(cache.algebraic_vars, fsallast, false) /
integrator.opts.abstol
invatol = inv(integrator.opts.abstol)
@.. atmp=ifelse(cache.algebraic_vars, fsallast, false) * invatol
integrator.EEst += integrator.opts.internalnorm(atmp, t)
end
end
Expand All @@ -244,7 +244,8 @@ end
@unpack c₃₂, d, tf, uf = cache

# Precalculations
γ = dt * d
dtγ = dt * d
neginvdtγ = -inv(dtγ)
dto2 = dt / 2
dto6 = dt / 6

Expand All @@ -258,22 +259,24 @@ end
# Time derivative
dT = calc_tderivative(integrator, cache)

W = calc_W(integrator, cache, γ, repeat_step)
W = calc_W(integrator, cache, dtγ, repeat_step, true)
if !issuccess_W(W)
integrator.EEst = 2
return nothing
end

k₁ = _reshape(W \ -_vec((integrator.fsalfirst + γ * dT)), axes(uprev))
k₁ = _reshape(W \ _vec((integrator.fsalfirst + dtγ * dT)), axes(uprev)) * neginvdtγ
integrator.stats.nsolve += 1
f₁ = f(uprev + dto2 * k₁, p, t + dto2)
tmp = @.. uprev + dto2 * k₁
f₁ = f(tmp, p, t + dto2)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)

if mass_matrix === I
k₂ = _reshape(W \ -_vec(f₁ - k₁), axes(uprev)) + k₁
k₂ = _reshape(W \ _vec(f₁ - k₁), axes(uprev))
else
k₂ = _reshape(W \ -_vec(f₁ - mass_matrix * k₁), axes(uprev)) + k₁
k₂ = _reshape(W \ _vec(f₁ - mass_matrix * k₁), axes(uprev))
end
k₂ = @.. k₂ * neginvdtγ + k₁
integrator.stats.nsolve += 1
u = uprev + dt * k₂

Expand All @@ -282,30 +285,28 @@ end
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)

if mass_matrix === I
k₃ = _reshape(
W \
-_vec((integrator.fsallast - c₃₂ * (k₂ - f₁) -
2 * (k₁ - integrator.fsalfirst) + dt * dT)),
axes(uprev))
linsolve_tmp = @.. (integrator.fsallast - c₃₂ * (k₂ - f₁) -
2 * (k₁ - integrator.fsalfirst) + dt * dT)
else
linsolve_tmp = integrator.fsallast - mass_matrix * (c₃₂ * k₂ + 2 * k₁) +
c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT
k₃ = _reshape(W \ -_vec(linsolve_tmp), axes(uprev))
linsolve_tmp = mass_matrix * (@.. c₃₂ * k₂ + 2 * k₁)
linsolve_tmp = @.. (integrator.fsallast - linsolve_tmp +
c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT)
end
k₃ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) * neginvdtγ
integrator.stats.nsolve += 1

if u isa Number
utilde = dto6 * f.mass_matrix[1, 1] * (k₁ - 2 * k₂ + k₃)
else
utilde = dto6 * f.mass_matrix * (k₁ - 2 * k₂ + k₃)
utilde = f.mass_matrix * (@.. dto6 * (k₁ - 2 * k₂ + k₃))
end
atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol,
integrator.opts.reltol, integrator.opts.internalnorm, t)
integrator.EEst = integrator.opts.internalnorm(atmp, t)

if mass_matrix !== I
atmp = @. ifelse(!integrator.differential_vars, integrator.fsallast, false) ./
integrator.opts.abstol
invatol = inv(integrator.opts.abstol)
atmp = @. ifelse(integrator.differential_vars, false, integrator.fsallast) * invatol
integrator.EEst += integrator.opts.internalnorm(atmp, t)
end
end
Expand All @@ -321,7 +322,8 @@ end
@unpack c₃₂, d, tf, uf = cache

# Precalculations
γ = dt * d
dtγ = dt * d
neginvdtγ = -inv(dtγ)
dto2 = dt / 2
dto6 = dt / 6

Expand All @@ -335,52 +337,52 @@ end
# Time derivative
dT = calc_tderivative(integrator, cache)

W = calc_W(integrator, cache, γ, repeat_step)
W = calc_W(integrator, cache, dtγ, repeat_step, true)
if !issuccess_W(W)
integrator.EEst = 2
return nothing
end

k₁ = _reshape(W \ -_vec((integrator.fsalfirst + γ * dT)), axes(uprev))
k₁ = _reshape(W \ -_vec((integrator.fsalfirst + dtγ * dT)), axes(uprev))/dtγ
integrator.stats.nsolve += 1
f₁ = f(uprev + dto2 * k₁, p, t + dto2)
tmp = @.. uprev + dto2 * k₁
f₁ = f(tmp, p, t + dto2)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)

if mass_matrix === I
k₂ = _reshape(W \ -_vec(f₁ - k₁), axes(uprev)) + k₁
k₂ = _reshape(W \ _vec(f₁ - k₁), axes(uprev))
else
linsolve_tmp = f₁ - mass_matrix * k₁
k₂ = _reshape(W \ -_vec(linsolve_tmp), axes(uprev)) + k₁
k₂ = _reshape(W \ _vec(linsolve_tmp), axes(uprev))
end
k₂ = @.. k₂ * neginvdtγ + k₁

integrator.stats.nsolve += 1
tmp = uprev + dt * k₂
tmp = @.. uprev + dt * k₂
integrator.fsallast = f(tmp, p, t + dt)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)

if mass_matrix === I
k₃ = _reshape(
W \
-_vec((integrator.fsallast - c₃₂ * (k₂ - f₁) -
2(k₁ - integrator.fsalfirst) + dt * dT)),
axes(uprev))
linsolve_tmp = @.. (integrator.fsallast - c₃₂ * (k₂ - f₁) -
2(k₁ - integrator.fsalfirst) + dt * dT)
else
linsolve_tmp = integrator.fsallast - mass_matrix * (c₃₂ * k₂ + 2k₁) + c₃₂ * f₁ +
2 * integrator.fsalfirst + dt * dT
k₃ = _reshape(W \ -_vec(linsolve_tmp), axes(uprev))
linsolve_tmp = mass_matrix * (@.. c₃₂ * k₂ + 2 * k₁)
linsolve_tmp = @.. (integrator.fsallast - linsolve_tmp +
c₃₂ * f₁ + 2 * integrator.fsalfirst + dt * dT)
end
k₃ = _reshape(W \ _vec(linsolve_tmp), axes(uprev)) * neginvdtγ
integrator.stats.nsolve += 1
u = uprev + dto6 * (k₁ + 4k₂ + k₃)
u = @.. uprev + dto6 * (k₁ + 4k₂ + k₃)

if integrator.opts.adaptive
utilde = dto6 * (k₁ - 2k₂ + k₃)
utilde = @.. dto6 * (k₁ - 2k₂ + k₃)
atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol,
integrator.opts.reltol, integrator.opts.internalnorm, t)
integrator.EEst = integrator.opts.internalnorm(atmp, t)

if mass_matrix !== I
atmp = @. ifelse(!integrator.differential_vars, integrator.fsallast, false) ./
integrator.opts.abstol
invatol = inv(integrator.opts.abstol)
atmp = ifelse(integrator.differential_vars, false, integrator.fsallast) .* invatol
integrator.EEst += integrator.opts.internalnorm(atmp, t)
end
end
Expand Down
Loading
Loading