From 5410dc627ea4690ad0734d7ec1eabb252e0d0e58 Mon Sep 17 00:00:00 2001 From: FHoltorf Date: Mon, 29 Jan 2024 23:53:23 -0500 Subject: [PATCH] adjust initialization for correct fsal updates --- src/integrators/on_the_fly_interpolation.jl | 2 +- src/integrators/projector_splitting.jl | 47 +++++++++------ .../rank_adaptive_unconventional.jl | 58 ++++++++++++++----- src/integrators/unconventional.jl | 30 ++++++---- 4 files changed, 93 insertions(+), 44 deletions(-) diff --git a/src/integrators/on_the_fly_interpolation.jl b/src/integrators/on_the_fly_interpolation.jl index df4cad2..7476af1 100644 --- a/src/integrators/on_the_fly_interpolation.jl +++ b/src/integrators/on_the_fly_interpolation.jl @@ -24,7 +24,7 @@ function select_idcs(F::SVDLikeRepresentation, selection_alg::IndexSelectionAlgo row_idcs = index_selection(F_svd.U, diag(F_svd.S), selection_alg) col_idcs = index_selection(F_svd.V, diag(F_svd.S), selection_alg) UF = F_svd.U[:, 1:length(row_idcs)] - VF = F_trunc.V[:, 1:length(col_idcs)] + VF = F_svd.V[:, 1:length(col_idcs)] return UF, VF, row_idcs, col_idcs end diff --git a/src/integrators/projector_splitting.jl b/src/integrators/projector_splitting.jl index 71c3a76..6290750 100644 --- a/src/integrators/projector_splitting.jl +++ b/src/integrators/projector_splitting.jl @@ -50,7 +50,6 @@ end function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = prob.tspan[1]) # allocate memory for frequently accessed arrays tspan = (t0,t0+dt) - US = u.U*u.S VS = u.V*u.S' QRK = qr(US) @@ -64,7 +63,9 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p K_rhs = alg.alg_params.K_rhs end KProblem = ODEProblem(K_rhs, US, tspan, u.V) - KIntegrator = init(KProblem, alg.alg_params.K_alg, save_everystep=false, alg.alg_params.K_kwargs...) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) if isnothing(alg.alg_params.S_rhs) S_rhs = function (S, (U,V), t) @@ -74,7 +75,9 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p S_rhs = alg.alg_params.S_rhs end SProblem = ODEProblem(S_rhs, QRK.R, tspan, (u.U, u.V)) - SIntegrator = init(SProblem, alg.alg_params.S_alg, save_everystep=false, alg.alg_params.S_kwargs...) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) if isnothing(alg.alg_params.L_rhs) L_rhs = function (VS, U, t) @@ -84,8 +87,10 @@ function alg_cache(prob::MatrixDEProblem, alg::ProjectorSplitting, u, dt; t0 = p L_rhs = alg.alg_params.L_rhs end LProblem = ODEProblem(L_rhs, VS, tspan, u.U) - LIntegrator = init(LProblem, alg.alg_params.L_alg, save_everystep=false, alg.alg_params.L_kwargs...) - + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + return ProjectorSplitting_Cache(US, VS, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, nothing) @@ -126,11 +131,13 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) - + step!(KIntegrator,0.01*dt,true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (dS, S, p, t) Π_S, U1, V1, params = p - -Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) + Π_S(dS, SVDLikeRepresentation(U1,S,V1), params, t) end else S_rhs = alg.alg_params.S_rhs @@ -143,7 +150,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, QRK.R, tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - + step!(SIntegrator,0.01*dt,true) + set_t!(SIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (dL, L, p, t) Π_L, U0, params = p @@ -158,10 +167,12 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_L = (Π_L, u.U, ()) LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) - + step!(LIntegrator,0.01*dt,true) + set_t!(LIntegrator, t0) + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), - Π, Π_K, Π_L, Π_S) + Π, Π_K, Π_L, Π_S) return ProjectorSplitting_Cache(US, VS, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, interpolation_cache) @@ -282,15 +293,16 @@ function primal_LT_deim_step!(u, cache, t, dt) set_u!(KIntegrator, US) step!(KIntegrator, dt, true) US .= KIntegrator.u - QRL = qr!(US) - u.U .= Matrix(QRL.Q) - mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights) - + QRK = qr!(US) + u.U .= Matrix(QRK.Q) + # S step - set_u!(SIntegrator, QRL.R) + mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights, -1.0, 0) + set_u!(SIntegrator, QRK.R) step!(SIntegrator, dt, true) # L step + mul!(Π_L.interpolator.weights, u.U', Π.interpolator.range.weights) mul!(VS,u.V,SIntegrator.u') set_u!(LIntegrator, VS) step!(LIntegrator, dt, true) @@ -325,14 +337,15 @@ function dual_LT_deim_step!(u, cache, t, dt) VS .= LIntegrator.u QRL = qr!(VS) u.V .= Matrix(QRL.Q) - mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights) - + # S step + mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights, -1, 0) set_u!(SIntegrator, Matrix(QRL.R')) step!(SIntegrator, dt, true) # K step mul!(US,u.U,SIntegrator.u) + mul!(Π_K.interpolator.parent.weights, u.V', Π.interpolator.corange.weights) set_u!(KIntegrator, US) step!(KIntegrator, dt, true) US .= KIntegrator.u diff --git a/src/integrators/rank_adaptive_unconventional.jl b/src/integrators/rank_adaptive_unconventional.jl index 9cbf25a..6539954 100644 --- a/src/integrators/rank_adaptive_unconventional.jl +++ b/src/integrators/rank_adaptive_unconventional.jl @@ -69,6 +69,9 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit end KProblem = ODEProblem(K_rhs, US, tspan, u.V) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (VS, U, t) return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U) @@ -76,9 +79,11 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit else L_rhs = alg.alg_params.L_rhs end - LProblem = ODEProblem(L_rhs, VS, tspan, u.U) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (S, (U,V), t) return Matrix(U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V) @@ -88,6 +93,8 @@ function alg_cache(prob::MatrixDEProblem, alg::RankAdaptiveUnconventionalAlgorit end SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, [Uhat, Vhat]) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem, @@ -142,6 +149,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) if isnothing(alg.alg_params.L_rhs) L_rhs = function (dL, L, p, t) @@ -157,6 +166,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_L = (Π_L, u.U, ()) LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) if isnothing(alg.alg_params.S_rhs) S_rhs = function (dS, S, p, t) @@ -174,6 +185,8 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), Π, Π_K, Π_L, Π_S) @@ -226,31 +239,44 @@ function alg_recache(cache::RankAdaptiveUnconventionalAlgorithm_Cache, alg::Rank # the following sequence can probably somehow be handled via dispatch in a simpler way if isnothing(KProblem) KIntegrator = MatrixDataIntegrator(Δy, US, I, u.V, 1) - elseif isnothing(interpolation_cache) - KProblem = remake(KProblem, u0 = US, p = u.V, tspan = (t, t+1.0)) - KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) else - KProblem = remake(KProblem, u0 = US, p = (interpolation_cache.Π_K, u.V, ()), tspan = (t, t+1.0)) - KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + if isnothing(interpolation_cache) + KProblem = remake(KProblem, u0 = US, p = u.V, tspan = (t, t+1.0)) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + else + KProblem = remake(KProblem, u0 = US, p = (interpolation_cache.Π_K, u.V, ()), tspan = (t, t+1.0)) + KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) + end + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) end + if isnothing(LProblem) LIntegrator = MatrixDataIntegrator(Δy', VS, I, u.U, 1) - elseif isnothing(interpolation_cache) - LProblem = remake(LProblem, u0 = VS, p = u.U, tspan = (t, t+1.0)) - LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) else - LProblem = remake(LProblem, u0 = VS, p = (interpolation_cache.Π_L, u.U, ()), tspan = (t, t+1.0)) - LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + if isnothing(interpolation_cache) + LProblem = remake(LProblem, u0 = VS, p = u.U, tspan = (t, t+1.0)) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + else + LProblem = remake(LProblem, u0 = VS, p = (interpolation_cache.Π_L, u.U, ()), tspan = (t, t+1.0)) + LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) + end + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) end if isnothing(SProblem) SIntegrator = MatrixDataIntegrator(Δy, zeros(2r,2r), Uhat, Vhat, 1) - elseif isnothing(interpolation_cache) - SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = [Uhat, Vhat], tspan = (t, t+1.0)) - SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) else - SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = (interpolation_cache.Π_S, U_hat, V_hat, ()), tspan = (t, t+1.0)) - SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + if isnothing(interpolation_cache) + SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = [Uhat, Vhat], tspan = (t, t+1.0)) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + else + SProblem = remake(SProblem, u0 = zeros(2*r,2*r), p = (interpolation_cache.Π_S, U_hat, V_hat, ()), tspan = (t, t+1.0)) + SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) + end + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) end return RankAdaptiveUnconventionalAlgorithm_Cache(US, Uhat, VS, Vhat, M, N, SIntegrator, SProblem, LIntegrator, LProblem, diff --git a/src/integrators/unconventional.jl b/src/integrators/unconventional.jl index 764be2e..cd9c687 100644 --- a/src/integrators/unconventional.jl +++ b/src/integrators/unconventional.jl @@ -80,7 +80,9 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t end KProblem = ODEProblem(K_rhs, US, tspan, u.V) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) - + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (VS, U, t) return Matrix(prob.f(TwoFactorRepresentation(U,VS),t)'*U) @@ -90,7 +92,9 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t end LProblem = ODEProblem(L_rhs, VS, tspan, u.U) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) - + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (S, (U,V), t) return Matrix(U'*prob.f(SVDLikeRepresentation(U,S,V),t)*V) @@ -100,17 +104,17 @@ function alg_cache(prob::MatrixDEProblem, alg::UnconventionalAlgorithm, u, dt; t end SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, (u.U, u.V)) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) + return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, nothing) end - function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, alg::UnconventionalAlgorithm, u, dt; t0 = prob.tspan[1]) where {fType <: ComponentFunction, uType, tType} - # allocate memory for frequently accessed arrays tspan = (t0,t0+dt) @unpack tol, rmin, rmax, @@ -146,7 +150,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_K = (Π_K, u.V, ()) KProblem = ODEProblem(K_rhs, US, tspan, p_K) KIntegrator = init(KProblem, alg.alg_params.K_alg; save_everystep=false, alg.alg_params.K_kwargs...) - + step!(KIntegrator, 0.01*dt, true) + set_t!(KIntegrator, t0) + if isnothing(alg.alg_params.L_rhs) L_rhs = function (dL, L, p, t) Π_L, U0, params = p @@ -161,7 +167,9 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, p_L = (Π_L, u.U, ()) LProblem = ODEProblem(L_rhs, VS, tspan, p_L) LIntegrator = init(LProblem, alg.alg_params.L_alg; save_everystep=false, alg.alg_params.L_kwargs...) - + step!(LIntegrator, 0.01*dt, true) + set_t!(LIntegrator, t0) + if isnothing(alg.alg_params.S_rhs) S_rhs = function (dS, S, p, t) Π_S, U1, V1, params = p @@ -178,9 +186,11 @@ function alg_cache(prob::MatrixDEProblem{fType, uType, tType}, SProblem = ODEProblem(S_rhs, M*u.S*N', tspan, p_S) SIntegrator = init(SProblem, alg.alg_params.S_alg; save_everystep=false, alg.alg_params.S_kwargs...) - + step!(SIntegrator, 0.01*dt, true) + set_t!(SIntegrator, t0) + interpolation_cache = OnTheFlyInterpolation_Cache(alg.alg_params.interpolation, deepcopy(u), - Π, Π_K, Π_L, Π_S) + Π, Π_K, Π_L, Π_S) return UnconventionalAlgorithm_Cache(US, VS, M, N, QRK, QRL, SIntegrator, LIntegrator, KIntegrator, nothing, nothing, nothing, nothing, interpolation_cache) @@ -282,7 +292,7 @@ function unconventional_deim_step!(u, cache, t, dt) mul!(N,Matrix(QRL.Q)',u.V) u.V .= Matrix(QRL.Q) u.U .= Matrix(QRK.Q) - + # update core interpolator mul!(Π_S.interpolator.range.weights, u.U', Π.interpolator.range.weights) mul!(Π_S.interpolator.corange.weights, u.V', Π.interpolator.corange.weights)