search

SKopecz / PositiveIntegrators.jl

A Julia library of positivity-preserving time integration methods
https://skopecz.github.io/PositiveIntegrators.jl/
MIT License
13 stars 4 forks source link

Careful handling of data types (Float32 & units) #108

Closed SKopecz closed 3 weeks ago

SKopecz commented 2 months ago
codecov-commenter commented 2 months ago

:warning: Please install the 'codecov app svg image' to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 97.72727% with 1 line in your changes missing coverage. Please review.

Project coverage is 98.26%. Comparing base (e9de366) to head (cba3acb).

Files with missing lines Patch % Lines
src/sspmprk.jl 92.30% 1 Missing :warning:

:exclamation: Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #108 +/- ## ========================================== + Coverage 98.24% 98.26% +0.02% ========================================== Files 7 7 Lines 1538 1556 +18 ========================================== + Hits 1511 1529 +18 Misses 27 27 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

SKopecz commented 2 months ago

The code

    linsolve.A = P3
    @show linsolve.A
    @show linsolve.b
    A_tmp = Matrix(linsolve.A)
    b_tmp = linsolve.b
    linres = solve!(linsolve)
    @show linres
    @show A_tmp\b_tmp
    integrator.stats.nsolve += 1

    σ .= linres
    @show ("6",σ)

provides

linsolve.A = sparse([1, 2, 1, 2, 3, 2, 3, 4, 3, 4], [1, 1, 2, 2, 2, 3, 3, 3, 4, 4], [1.0, -0.0, -0.0001306202464031375, 1.000130620246403, -0.0, -6.4983581002939935e59, 6.4983581002939935e59, -0.0, -3.400386856217244e261, 3.400386856217244e261], 4, 4)
linsolve.b = [2.2436232870583095, 4.756376712918744, 3.6441102462903125e-43, 0.0]
linres = [2.244244485015507, 4.755755514961547, -5.861544385515315e-76, 0.0]
A_tmp \ b_tmp = [2.244244485015507, 4.755755514961547, 5.60773997069421e-103, 0.0]
("6", σ) = ("6", [2.244244485015507, 4.755755514961547, -5.861544385515315e-76, 0.0])

So the problem is indeed a linear solve which introduces a negative element in the solution.

SKopecz commented 2 months ago

@ranocha How shall we deal with this?

It is clear on the analytical level that the solutions of the linear system should be nonnegative, since the matrix A is an M-matrix and the right hand side b is nonnegative. But I wonder if this is actually guaranteed on the numerical level?

ranocha commented 2 months ago

I don't think there is a simple way to guarantee this at the numerical level. Even small floating point errors can change a value that is zero to something slightly negative. In your example above, -5.861544385515315e-76 is really tiny but negative compared to the order unity values.

ranocha commented 2 months ago

We could try to play with small_constant for this specific test case or change the problem setup slightly to avoid the error in CI here - or continue ignoring it. However, this does not fix the underlying problem.

SKopecz commented 2 months ago

Setting small_constant = 1e-50 worked indeed. This is the same value which is used in SSPMPRK43().

SKopecz commented 1 month ago

The test below shows that for almost all schemes the constants have the same type as the algorithm input parameters. So all computations in perform_step! should be in single precision, whenever Float32 arguments are passed to an algorithm.

The only exception is SSPMPRK43() which takes no inputs and has hard codedFloat64 constants. So we only need to take special care of this algorithm.

@testset "Float32" begin
    algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0),
            MPRK43I(0.5f0, 0.75f0), MPRK43II(0.5f0), MPRK43II(2.0f0 / 3.0f0),
            SSPMPRK22(0.5f0, 1.0f0), SSPMPRK43())
    @testset "$i" for (i, alg) in enumerate(algs)
        @test eltype(PositiveIntegrators.get_constant_parameters(alg)) == Float32
    end

    algs = (MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5),
            MPRK43I(0.5, 0.75), MPRK43II(0.5), MPRK43II(2.0 / 3.0),
            SSPMPRK22(0.5, 1.0), SSPMPRK43())
    @testset "$i" for (i, alg) in enumerate(algs)
        @test eltype(PositiveIntegrators.get_constant_parameters(alg)) == Float64
    end
end
Test Summary: | Pass  Fail  Total  Time
Float32       |   17     1     18  0.1s
  1           |    1            1  0.0s
  2           |    1            1  0.0s
  3           |    1            1  0.0s
  4           |    1            1  0.0s
  5           |    1            1  0.0s
  6           |    1            1  0.0s
  7           |    1            1  0.0s
  8           |    1            1  0.0s
  9           |          1      1  0.0s
  1           |    1            1  0.0s
  2           |    1            1  0.0s
  3           |    1            1  0.0s
  4           |    1            1  0.1s
  5           |    1            1  0.0s
  6           |    1            1  0.0s
  7           |    1            1  0.0s
  8           |    1            1  0.0s
  9           |    1            1  0.0s
ranocha commented 1 month ago

The test below shows that for almost all schemes the constants have the same type as the algorithm input parameters. So all computations in perform_step! should be in single precision, whenever Float32 arguments are passed to an algorithm.

The only exception is SSPMPRK43() which takes no inputs and has hard codedFloat64 constants. So we only need to take special care of this algorithm.

We call get_constant_parameters(alg) in alg_cache, where we have access to uBottomEltypeNoUnits and tTypeNoUnits. Can we pass these to get_constant_parameters like in get_constant_parameters(alg, uBottomEltypeNoUnits, tTypeNoUnits) and then use them to convert the coefficients appropriately, e.g., 

function get_constant_parameters(alg::MPRK22, ::Type{Tu}, ::Type{Tt}) where {Tu, Tt}
    # ...

    a21 = convert(Tu, alg.alpha)
    b2 = convert(T0, 1 / (2 * a21))
    b1 = convert(Tu, 1 - b2)

instead of

function get_constant_parameters(alg::MPRK22)
    # ...

    a21 = alg.alpha
    b2 = 1 / (2 * a21)
    b1 = 1 - b2

for MPRK22 and similar for the other algorithms? Coefficients that would only multiply times would use Tt and the other ones Tu.

SKopecz commented 1 month ago

What's the thinking behind this? Why would we want Float32 times and Float64 approximations u? Or are there other data types for which this would be of interest?

ranocha commented 1 month ago

It's likely not something that we will use here, but OrdinaryDiffEq.jl allows using types with units or complex numbers for u

ranocha commented 1 month ago

I think we can simplify it here and just pass uBottomEltypeNoUnits when setting up the parameters during construction of the cache

SKopecz commented 1 month ago

Complex numbers in the world of PDS is an unsolved problem I guess. But units might be of interest. I'll see if I can make it work.

What exactly is the difference between uEltypeNoUnits and uBottomEltypeNoUnits?

I assume small_constant needs to be of uType then.

ranocha commented 1 month ago

The element type of u can be complex or dual numbers. In that case, the bottom type is the underlying real type

ranocha commented 1 month ago

Amd yes, the small constant should be of the element type of u

SKopecz commented 3 weeks ago

This PR contains much more than originally planned.

SKopecz commented 3 weeks ago

I came across https://discourse.julialang.org/t/unitful-jl-and-differentialequations-jl-compatibility/68594, where it is discussed that Unitful.jl cannot deal with linear algebra. As all MPRK schemes require linear solves, these schemes cannot deal with units from Unitful.jl. The same is true for implicit solvers from OrdinaryDiffEq.

The post also mentions DynamicQuantities.jl or UnitfulLinearAlgebra, but I think we should focus on other things. So this PR could be merged.

ranocha commented 3 weeks ago

The Downgrade tests fail. Do you have an idea why?

SKopecz commented 3 weeks ago

The Downgrade tests fail. Do you have an idea why?

Is there a way to reproduce these locally?

ranocha commented 3 weeks ago

You can create Project.toml files with the versions of the direct dependencies listed in https://github.com/SKopecz/PositiveIntegrators.jl/actions/runs/10580596414/job/29315860050?pr=108#step:6:11

SKopecz commented 3 weeks ago

The Downgrade tests fail. Do you have an idea why?

The downgrade tests didn't accept isapprox(sol1, sol2) for solutions with units. Now we use isapprox(ustrip.(sol1.u),ustrip.(sol2.u)).