Careful handling of data types (Float32 & units)

[SKopecz]

• Display variables during perform_step

[codecov-commenter]

⚠️ Please install the 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 ⚠️

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

Additional details and impacted files

@@            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              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[SKopecz]

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]

@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]

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]

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]

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

[SKopecz]

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]

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]

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]

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

[ranocha]

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

[SKopecz]

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]

Thanks! Please find below some suggestions

[ranocha]

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

[SKopecz]

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

Is there a way to reproduce these locally?

[ranocha]

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]

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)).

[ranocha]

Thanks!