More flexibility for internal types

[termi-official]

This PR fixes two weakly related things.

First, it allows to actually pass the matrix type into the sparsity pattern constructor. Currently only SparseMatrixCSC with element type of Float64 was allowed. Now we can also use different element types. Other matrix types are not yet implemented and up to future PRs (also related to FerriteDistributed.jl).

The second thing is that at some code points we fixed the data types to Float64. This constraint has also been removed, allowing e.g. the usage of Float32. This change also opens the possibility to implement an example of a complex-valued problem (e.g. Schrödinger problem).

This PR precedes #766 , because it can be useful to use lower precision on the GPU, e.g. to increase throughput and because not every GPU has dedicated Float64 units.

Closes #195 .

[KnutAM]

Nice!
Just leaving the note here, not sure if already on the list to change, that many of the shape_value(ip, \xi) implementations hard-code Float64, but should be easy to change (might also not matter as it in many cases would convert to a container type of lower precision), e.g.

Ferrite.jl/src/interpolations.jl

Line 535 in 383c0c5

function shape_value(ip::Lagrange{RefQuadrilateral, 1}, ξ::Vec{2}, i::Int)

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (324bd4a) 93.28% compared to head (4843a07) 93.28%.
Report is 3 commits behind head on master.

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

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #802      +/-   ##
==========================================
- Coverage   93.28%   93.28%   -0.01%     
==========================================
  Files          36       36              
  Lines        5228     5224       -4     
==========================================
- Hits         4877     4873       -4     
  Misses        351      351              

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

[termi-official]

Yea, this is why I have left it as a draft. I am not sure how to resolve this one properly, because there are several possibilities. First, we could make the value type part of e.g. the Lagrange struct parametrization. On the other hand, we also might be able to have extra dispatches on e.g. shape_value. Not entirely sure. But the body should be definitely in the correct precision!

[KnutAM]

I think in most cases it should be as simple as

function shape_value(ip::Lagrange{RefQuadrilateral, 1}, ξ::Vec{2}, i::Int)
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return (1 - ξ_x) * (1 - ξ_y) * 0.25
    i == 2 && return (1 + ξ_x) * (1 - ξ_y) * 0.25
    i == 3 && return (1 + ξ_x) * (1 + ξ_y) * 0.25
    i == 4 && return (1 - ξ_x) * (1 + ξ_y) * 0.25
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

to

function shape_value(ip::Lagrange{RefQuadrilateral, 1}, ξ::Vec{2,T}, i::Int) where T
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 1 && return (oneunit(T) - ξ_x) * (oneunit(T) - ξ_y) / 4
    i == 2 && return (oneunit(T) + ξ_x) * (oneunit(T) - ξ_y) / 4
    i == 3 && return (oneunit(T) + ξ_x) * (oneunit(T) + ξ_y) / 4
    i == 4 && return (oneunit(T) - ξ_x) * (oneunit(T) + ξ_y) / 4
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

(I suppose this is what you meant with the body, but left it here anyways)

If supplied a quadrature rule with the correct precision/type of the xi-coordinates, I suppose it should work? But there are probably more entry points in some cases, such as the reference coordinates when using BCValues for Dirichlet boundary conditions.

[termi-official]

Okay, this took longer than it should have taken. :D Interestingly 1 seems to work fine.

[termi-official]

Following the suggestion by Knut in #766 (comment) I also included a change from ArgumentError to BoundsError for interpolations (which might fix the issue in the linked PR).

[KnutAM]

Interestingly 1 seems to work fine.

Using Int would work as long as there are no units, but this would throw if someone tries to supply coordinate with units - but this will break a lot of other assumptions in the code base too, so perhaps not worth the extra hassle :)

[termi-official]

Interestingly 1 seems to work fine.

Using Int would work as long as there are no units, but this would throw if someone tries to supply coordinate with units - but this will break a lot of other assumptions in the code base too, so perhaps not worth the extra hassle :)

I think this is still interesting to explore. :)

[KnutAM]

Nice work! Quite a lot of detailed changes!

I've tried to go through carefully but especially in the ConstraintHandler, there seem to be several possibilities for Float64 to sneak in.
Not sure if it is necessary to prevent this as long as (a) we have type-stable code and (2) the output types are the correct types. I guess further improvements can be done in the GPU PR if issues arise with internal Float64 then?

In addition to the JET testing, it would be nice to have a few @btime comparisons with master for usages of the ConstraintHandler with Float64.
For example

• close!
• update!
• apply!
• Creation of periodic boundary conditions
  Where the 3 first are done for cases including affine constraints to avoid any regressions here. This could also, in addition, be tested with JET, if possible.

[termi-official]

Thanks for the detailed review!

I've tried to go through carefully but especially in the ConstraintHandler, there seem to be several possibilities for Float64 to sneak in. Not sure if it is necessary to prevent this as long as (a) we have type-stable code and (2) the output types are the correct types. I guess further improvements can be done in the GPU PR if issues arise with internal Float64 then?

I think I would like to fix such issues in this PR.

In addition to the JET testing, it would be nice to have a few @btime comparisons with master for usages of the ConstraintHandler with Float64. For example

* `close!`

* `update!`

* `apply!`

* Creation of periodic boundary conditions
  Where the 3 first are done for cases including affine constraints to avoid any regressions here. This could also, in addition, be tested with JET, if possible.

using BenchmarkTools, Ferrite

grid = generate_grid(Quadrilateral, (1000, 1000));
ip = Lagrange{RefQuadrilateral, 1}()
qr = QuadratureRule{RefQuadrilateral}(2)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh);

@btime ConstraintHandler(dh)

∂Ω = union(
    getfaceset(grid, "left"),
    getfaceset(grid, "right"),
    getfaceset(grid, "top"),
    getfaceset(grid, "bottom"),
);
dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
@btime add!(ch, $dbc) setup=(ch = ConstraintHandler($dh)) evals=1;

function setup_ch(dh, dbc)
    ch = ConstraintHandler(dh)
    add!(ch, dbc)
    return ch
end
@btime close!(ch) setup=(ch = setup_ch($dh, $dbc)) evals=1;

function setup_ch(dh, dbc)
    ch = ConstraintHandler(dh)
    add!(ch, dbc)
    close!(ch)
    return ch
end
@btime update!(ch) setup=(ch = setup_ch($dh, $dbc)) evals=1;

shows only a difference in the ctor.

[termi-official]

So I guess #818 precedes this. :D