|
88 | 88 | sol2 = solve(prob2, Tsit5()) |
89 | 89 | @test 0.0:0.07:1.0 ⊆ sol2.t |
90 | 90 | end |
| 91 | + |
| 92 | +@testset "Tstop Robustness Tests (StaticArrays vs Arrays)" begin |
| 93 | + # Tests for issue #2752: tstop overshoot errors with StaticArrays |
| 94 | + |
| 95 | + @testset "StaticArrays vs Arrays with extreme precision" begin |
| 96 | + # Test the specific case that was failing: extreme precision + StaticArrays |
| 97 | + function precise_dynamics(u, p, t) |
| 98 | + x = @view u[1:2] |
| 99 | + v = @view u[3:4] |
| 100 | + |
| 101 | + # Electromagnetic-like dynamics |
| 102 | + dv = -0.01 * x + 1e-6 * sin(100*t) * SVector{2}(1, 1) |
| 103 | + |
| 104 | + return SVector{4}(v[1], v[2], dv[1], dv[2]) |
| 105 | + end |
| 106 | + |
| 107 | + function precise_dynamics_array!(du, u, p, t) |
| 108 | + x = @view u[1:2] |
| 109 | + v = @view u[3:4] |
| 110 | + |
| 111 | + dv = -0.01 * x + 1e-6 * sin(100*t) * [1, 1] |
| 112 | + du[1] = v[1] |
| 113 | + du[2] = v[2] |
| 114 | + du[3] = dv[1] |
| 115 | + du[4] = dv[2] |
| 116 | + end |
| 117 | + |
| 118 | + # Initial conditions |
| 119 | + u0_static = SVector{4}(1.0, -0.5, 0.01, 0.01) |
| 120 | + u0_array = [1.0, -0.5, 0.01, 0.01] |
| 121 | + tspan = (0.0, 2.0) |
| 122 | + tstops = [0.5, 1.0, 1.5] |
| 123 | + |
| 124 | + # Test with extreme tolerances that originally caused issues |
| 125 | + prob_static = ODEProblem(precise_dynamics, u0_static, tspan) |
| 126 | + sol_static = solve(prob_static, Vern9(); reltol=1e-12, abstol=1e-15, |
| 127 | + tstops=tstops, save_everystep=false) |
| 128 | + @test sol_static.retcode == :Success |
| 129 | + for tstop in tstops |
| 130 | + @test tstop ∈ sol_static.t |
| 131 | + end |
| 132 | + |
| 133 | + prob_array = ODEProblem(precise_dynamics_array!, u0_array, tspan) |
| 134 | + sol_array = solve(prob_array, Vern9(); reltol=1e-12, abstol=1e-15, |
| 135 | + tstops=tstops, save_everystep=false) |
| 136 | + @test sol_array.retcode == :Success |
| 137 | + for tstop in tstops |
| 138 | + @test tstop ∈ sol_array.t |
| 139 | + end |
| 140 | + |
| 141 | + # Solutions should be very close despite different array types |
| 142 | + @test isapprox(sol_static(2.0), sol_array(2.0), rtol=1e-10) |
| 143 | + end |
| 144 | + |
| 145 | + @testset "Duplicate tstops handling" begin |
| 146 | + function simple_ode(u, p, t) |
| 147 | + [0.1 * u[1]] |
| 148 | + end |
| 149 | + |
| 150 | + u0 = SVector{1}(1.0) |
| 151 | + tspan = (0.0, 2.0) |
| 152 | + |
| 153 | + # Test multiple identical tstops - should all be processed |
| 154 | + duplicate_tstops = [0.5, 0.5, 0.5, 1.0, 1.0] |
| 155 | + |
| 156 | + prob = ODEProblem(simple_ode, u0, tspan) |
| 157 | + sol = solve(prob, Vern9(); tstops=duplicate_tstops) |
| 158 | + |
| 159 | + @test sol.retcode == :Success |
| 160 | + |
| 161 | + # Count how many times each tstop appears in solution |
| 162 | + count_05 = count(t -> abs(t - 0.5) < 1e-12, sol.t) |
| 163 | + count_10 = count(t -> abs(t - 1.0) < 1e-12, sol.t) |
| 164 | + |
| 165 | + # Should handle all duplicate tstops (though may not save all due to deduplication) |
| 166 | + @test count_05 >= 1 # At least one 0.5 |
| 167 | + @test count_10 >= 1 # At least one 1.0 |
| 168 | + |
| 169 | + # Test with StaticArrays too |
| 170 | + prob_static = ODEProblem(simple_ode, u0, tspan) |
| 171 | + sol_static = solve(prob_static, Vern9(); tstops=duplicate_tstops) |
| 172 | + @test sol_static.retcode == :Success |
| 173 | + end |
| 174 | + |
| 175 | + @testset "PresetTimeCallback with identical times" begin |
| 176 | + # Test PresetTimeCallback scenarios where callbacks are set at same times as tstops |
| 177 | + |
| 178 | + event_times = Float64[] |
| 179 | + callback_times = Float64[] |
| 180 | + |
| 181 | + function affect_preset!(integrator) |
| 182 | + push!(callback_times, integrator.t) |
| 183 | + integrator.u[1] += 0.1 # Small modification |
| 184 | + end |
| 185 | + |
| 186 | + function simple_growth(u, p, t) |
| 187 | + [0.1 * u[1]] |
| 188 | + end |
| 189 | + |
| 190 | + u0 = SVector{1}(1.0) |
| 191 | + tspan = (0.0, 3.0) |
| 192 | + |
| 193 | + # Define times where both tstops and callbacks should trigger |
| 194 | + critical_times = [0.5, 1.0, 1.5, 2.0, 2.5] |
| 195 | + |
| 196 | + # Create PresetTimeCallback at the same times as tstops |
| 197 | + preset_cb = PresetTimeCallback(critical_times, affect_preset!) |
| 198 | + |
| 199 | + prob = ODEProblem(simple_growth, u0, tspan) |
| 200 | + sol = solve(prob, Vern9(); tstops=critical_times, callback=preset_cb, |
| 201 | + reltol=1e-10, abstol=1e-12) |
| 202 | + |
| 203 | + @test sol.retcode == :Success |
| 204 | + |
| 205 | + # Verify all tstops were hit |
| 206 | + for time in critical_times |
| 207 | + @test any(abs.(sol.t .- time) .< 1e-10) |
| 208 | + end |
| 209 | + |
| 210 | + # Verify all callbacks were triggered |
| 211 | + @test length(callback_times) == length(critical_times) |
| 212 | + for time in critical_times |
| 213 | + @test any(abs.(callback_times .- time) .< 1e-10) |
| 214 | + end |
| 215 | + |
| 216 | + # Test the same with regular arrays |
| 217 | + u0_array = [1.0] |
| 218 | + callback_times_array = Float64[] |
| 219 | + |
| 220 | + function affect_preset_array!(integrator) |
| 221 | + push!(callback_times_array, integrator.t) |
| 222 | + integrator.u[1] += 0.1 |
| 223 | + end |
| 224 | + |
| 225 | + function simple_growth_array!(du, u, p, t) |
| 226 | + du[1] = 0.1 * u[1] |
| 227 | + end |
| 228 | + |
| 229 | + preset_cb_array = PresetTimeCallback(critical_times, affect_preset_array!) |
| 230 | + |
| 231 | + prob_array = ODEProblem(simple_growth_array!, u0_array, tspan) |
| 232 | + sol_array = solve(prob_array, Vern9(); tstops=critical_times, callback=preset_cb_array, |
| 233 | + reltol=1e-10, abstol=1e-12) |
| 234 | + |
| 235 | + @test sol_array.retcode == :Success |
| 236 | + @test length(callback_times_array) == length(critical_times) |
| 237 | + |
| 238 | + # Both should have triggered all events |
| 239 | + @test length(callback_times) == length(callback_times_array) == length(critical_times) |
| 240 | + end |
| 241 | + |
| 242 | + @testset "Tiny tstop step handling" begin |
| 243 | + # Test cases where tstop is very close to current time (dt < eps(t)) |
| 244 | + function test_ode(u, p, t) |
| 245 | + [0.01 * u[1]] |
| 246 | + end |
| 247 | + |
| 248 | + u0 = SVector{1}(1.0) |
| 249 | + tspan = (0.0, 1.0) |
| 250 | + |
| 251 | + # Create tstop very close to start time (would cause tiny dt) |
| 252 | + tiny_tstops = [1e-15, 1e-14, 1e-13] |
| 253 | + |
| 254 | + for tiny_tstop in tiny_tstops |
| 255 | + prob = ODEProblem(test_ode, u0, tspan) |
| 256 | + sol = solve(prob, Vern9(); tstops=[tiny_tstop], save_everystep=false) |
| 257 | + |
| 258 | + @test sol.retcode == :Success |
| 259 | + @test any(abs.(sol.t .- tiny_tstop) .< 1e-14) # Should handle tiny tstop correctly |
| 260 | + end |
| 261 | + end |
| 262 | + |
| 263 | + @testset "Multiple close tstops with StaticArrays" begin |
| 264 | + # Test with multiple tstops that are very close together - stress test the flag logic |
| 265 | + function oscillator(u, p, t) |
| 266 | + SVector{2}(u[2], -u[1]) # Simple harmonic oscillator |
| 267 | + end |
| 268 | + |
| 269 | + u0 = SVector{2}(1.0, 0.0) |
| 270 | + tspan = (0.0, 4.0) |
| 271 | + |
| 272 | + # Multiple tstops close together (within floating-point precision range) |
| 273 | + close_tstops = [1.0, 1.0 + 1e-14, 1.0 + 2e-14, 1.0 + 5e-14, |
| 274 | + 2.0, 2.0 + 1e-15, 2.0 + 1e-14, |
| 275 | + 3.0, 3.0 + 1e-13] |
| 276 | + |
| 277 | + prob = ODEProblem(oscillator, u0, tspan) |
| 278 | + sol = solve(prob, Vern9(); tstops=close_tstops, reltol=1e-12, abstol=1e-15) |
| 279 | + |
| 280 | + @test sol.retcode == :Success |
| 281 | + |
| 282 | + # Should handle all close tstops without error |
| 283 | + # (Some might be deduplicated, but no errors should occur) |
| 284 | + unique_times = [1.0, 2.0, 3.0] |
| 285 | + for time in unique_times |
| 286 | + @test any(abs.(sol.t .- time) .< 1e-10) # At least hit the main times |
| 287 | + end |
| 288 | + end |
| 289 | + |
| 290 | + @testset "Backward integration with tstop flags" begin |
| 291 | + # Test that the fix works for backward time integration |
| 292 | + function decay_ode(u, p, t) |
| 293 | + [-0.1 * u[1]] |
| 294 | + end |
| 295 | + |
| 296 | + u0 = SVector{1}(1.0) |
| 297 | + tspan = (2.0, 0.0) # Backward integration |
| 298 | + tstops = [1.5, 1.0, 0.5] |
| 299 | + |
| 300 | + prob = ODEProblem(decay_ode, u0, tspan) |
| 301 | + sol = solve(prob, Vern9(); tstops=tstops, reltol=1e-12, abstol=1e-15) |
| 302 | + |
| 303 | + @test sol.retcode == :Success |
| 304 | + for tstop in tstops |
| 305 | + @test tstop ∈ sol.t |
| 306 | + end |
| 307 | + end |
| 308 | + |
| 309 | + @testset "Continuous callbacks during tstop steps" begin |
| 310 | + # Test that continuous callbacks work properly with tstop flag mechanism |
| 311 | + |
| 312 | + crossing_times = Float64[] |
| 313 | + |
| 314 | + function affect_continuous!(integrator) |
| 315 | + push!(crossing_times, integrator.t) |
| 316 | + end |
| 317 | + |
| 318 | + function condition_continuous(u, t, integrator) |
| 319 | + u[1] - 0.5 # Crosses when u[1] = 0.5 |
| 320 | + end |
| 321 | + |
| 322 | + function exponential_growth(u, p, t) |
| 323 | + [0.2 * u[1]] # Exponential growth |
| 324 | + end |
| 325 | + |
| 326 | + u0 = SVector{1}(0.1) # Start below 0.5 |
| 327 | + tspan = (0.0, 10.0) |
| 328 | + tstops = [2.0, 4.0, 6.0, 8.0] # Regular tstops |
| 329 | + |
| 330 | + continuous_cb = ContinuousCallback(condition_continuous, affect_continuous!) |
| 331 | + |
| 332 | + prob = ODEProblem(exponential_growth, u0, tspan) |
| 333 | + sol = solve(prob, Vern9(); tstops=tstops, callback=continuous_cb, |
| 334 | + reltol=1e-10, abstol=1e-12) |
| 335 | + |
| 336 | + @test sol.retcode == :Success |
| 337 | + |
| 338 | + # Should hit all tstops |
| 339 | + for tstop in tstops |
| 340 | + @test tstop ∈ sol.t |
| 341 | + end |
| 342 | + |
| 343 | + # Should also detect continuous callback crossings |
| 344 | + @test length(crossing_times) > 0 # At least one crossing detected |
| 345 | + |
| 346 | + # Verify crossings are at correct value |
| 347 | + for crossing_time in crossing_times |
| 348 | + u_at_crossing = sol(crossing_time) |
| 349 | + @test abs(u_at_crossing[1] - 0.5) < 1e-8 # Should be very close to 0.5 |
| 350 | + end |
| 351 | + end |
| 352 | + |
| 353 | +end |
0 commit comments