adjust initialization for correct fsal updates

FHoltorf · Jan 30, 2024 · 5410dc6 · 5410dc6
1 parent e0eba3d
commit 5410dc6
Show file tree

Hide file tree

Showing 4 changed files with 93 additions and 44 deletions.
diff --git 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
@@ -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
@@ -69,16 +69,21 @@ 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)
                 end
     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
@@ -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)