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

update #5

Merged
merged 4 commits into from
Jul 25, 2024
Merged

update #5

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7592944
don't broadcast the left side for out of place
oscardssmith Jul 24, 2024
367c34b
don't break error detection
oscardssmith Jul 24, 2024
17bd407
Merge pull request #2293 from oscardssmith/os/fix-oop-BDF
oscardssmith Jul 25, 2024
2020765
Merge pull request #4 from Shreyas-Ekanathan/master
Shreyas-Ekanathan Jul 25, 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
26 changes: 11 additions & 15 deletions lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,12 +47,11 @@ end
mass_matrix = f.mass_matrix

if mass_matrix === I
nlsolver.tmp = @.. broadcast=false ((dtₙ * β₁) * fₙ₋₁ +
(α₁ * uₙ₋₁ + α₂ * uₙ₋₂))/(dtₙ *
β₀)
nlsolver.tmp = @.. ((dtₙ * β₁) * fₙ₋₁ +
(α₁ * uₙ₋₁ + α₂ * uₙ₋₂))/(dtₙ * β₀)
else
_tmp = mass_matrix * @.. broadcast=false (α₁ * uₙ₋₁+α₂ * uₙ₋₂)
nlsolver.tmp = @.. broadcast=false ((dtₙ * β₁) * fₙ₋₁ + _tmp)/(dtₙ * β₀)
_tmp = mass_matrix * @.. (α₁ * uₙ₋₁+α₂ * uₙ₋₂)
nlsolver.tmp = @.. ((dtₙ * β₁) * fₙ₋₁ + _tmp)/(dtₙ * β₀)
end
nlsolver.γ = β₀
nlsolver.α = α₀
Expand Down Expand Up @@ -772,17 +771,17 @@ function perform_step!(integrator, cache::QNDFConstantCache{max_order},
u₀ = reshape(sum(D[:, 1:k], dims = 2) .+ uprev, size(u))
ϕ = zero(u)
for i in 1:k
@.. broadcast=false ϕ+=γₖ[i] * D[:, i]
ϕ = @.. ϕ + γₖ[i] * D[:, i]
end
end
markfirststage!(nlsolver)
nlsolver.z = u₀
mass_matrix = f.mass_matrix

if mass_matrix === I
nlsolver.tmp = @.. broadcast=false (u₀ / β₀ - ϕ)/dt
nlsolver.tmp = @.. (u₀ / β₀ - ϕ)/dt
else
nlsolver.tmp = mass_matrix * @.. broadcast=false (u₀ / β₀ - ϕ)/dt
nlsolver.tmp = mass_matrix * @.. (u₀ / β₀ - ϕ)/dt
end

nlsolver.γ = β₀
Expand Down Expand Up @@ -1110,18 +1109,16 @@ function perform_step!(integrator, cache::FBDFConstantCache{max_order},
end
else
for i in 1:(k - 1)
@.. broadcast=false @views u_corrector[:, i] = $calc_Lagrange_interp(
@.. @views u_corrector[:, i] = $calc_Lagrange_interp(
k, weights,
equi_ts[i],
ts,
u_history,
u_corrector[:,
i])
u_corrector[:, i])
end
tmp = -uprev * bdf_coeffs[k, 2]
vc = _vec(tmp)
for i in 1:(k - 1)
@.. broadcast=false @views vc -= u_corrector[:, i] * bdf_coeffs[k, i + 2]
@views tmp = @.. tmp - u_corrector[:, i] * bdf_coeffs[k, i + 2]
end
end

Expand Down Expand Up @@ -1177,9 +1174,8 @@ function perform_step!(integrator, cache::FBDFConstantCache{max_order},
end
terk *= abs(dt^(k))
else
vc = _vec(terk)
for i in 2:(k + 1)
@.. broadcast=false @views vc += fd_weights[i, k + 1] * u_history[:, i - 1]
@views terk = @.. terk + fd_weights[i, k + 1] * u_history[:, i - 1]
end
terk *= abs(dt^(k))
end
Expand Down
Loading