Helper script for convergence rates for Poisson problems

[termi-official]

I often use such scripts to test if my implementation converges correctly. We likely want to keep full convergence testing for all elements out of default CI since this can take an awful amount of time. Can be seen as integration level testing towards #558 .

[termi-official]

Alternatively I can also port over this script https://github.com/Ferrite-FEM/FerriteViz.jl/blob/master/docs/src/ferrite-examples/heat-equation.jl

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (e9ae940) 92.77% compared to head (6ce4346) 92.83%.
Report is 5 commits behind head on master.

❗ Current head 6ce4346 differs from pull request most recent head 86d98ff. Consider uploading reports for the commit 86d98ff to get more accurate results

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

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #640      +/-   ##
==========================================
+ Coverage   92.77%   92.83%   +0.06%     
==========================================
  Files          33       33              
  Lines        4952     4953       +1     
==========================================
+ Hits         4594     4598       +4     
+ Misses        358      355       -3     

Files	Coverage Δ
src/interpolations.jl	97.01% <100.00%> (+<0.01%)	⬆️

... and 1 file with indirect coverage changes

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

[KnutAM]

I find such an integration test to be a good idea. In order to catch a few more potential errors, perhaps it would be beneficial to have a case with fixed known solution which does not depend on calculating the rhs (I think e.g. errors in the numerical integration weights would not be caught by this test). Could perhaps be included in a github action to be run before each new release?

Did you intend to cast it into something like this?

Copied your code and played around a bit

The main point is the loop at the bottom...

using Ferrite, SparseArrays, ForwardDiff, Test
import LinearAlgebra: diag

struct TestCase{CT,O}
    qe_order::Int
    num_elems::Int
    point_tol::Float64
    L2_tol::Float64
    L∞_tol::Float64
end
function TestCase(;element, ip_order, qe_order, num_elems, point_tol=1e-2, L2_tol=1e-2, L∞_tol=1e-2)
    return TestCase{element, ip_order}(qe_order, num_elems, point_tol, L2_tol, L∞_tol)
end

# Solution parameters
analytical_solution(x) = prod(cos, x*π/2)
# Get the RHS analytically
analytical_rhs(x) = -sum(diag(ForwardDiff.hessian(analytical_solution,x)))

function setup_problem(case::TestCase{element,ip_order}) where {element, ip_order}
    qe_order = case.qe_order
    N = case.num_elems

    # Construct Ferrite stuff
    ip_geo = Ferrite.default_interpolation(element)
    dim = Ferrite.getdim(ip_geo)
    grid = generate_grid(element, ntuple(x->N, dim));
    refshape = Ferrite.getrefshape(ip_geo)
    ip = Lagrange{dim, refshape, ip_order}()
    qr = QuadratureRule{dim, refshape}(qe_order)
    cellvalues = CellScalarValues(qr, ip, ip_geo);
    dh = DofHandler(grid)
    add!(dh, :u, 1, ip)
    close!(dh);
    ch = ConstraintHandler(dh);
    ∂Ω = union(
        values(getfacesets(grid))...
    );
    dbc = Dirichlet(:u, ∂Ω, (x, t) -> analytical_solution(x))
    add!(ch, dbc);
    close!(ch)
    return dh, ch, cellvalues 
end

# Standard assembly copy pasta for Poisson problem
function assemble_element!(Ke::Matrix, fe::Vector, cellvalues::CellScalarValues, coords)
    n_basefuncs = getnbasefunctions(cellvalues)
    ## Reset to 0
    fill!(Ke, 0)
    fill!(fe, 0)
    ## Loop over quadrature points
    for q_point in 1:getnquadpoints(cellvalues)
        ## Get the quadrature weight
        dΩ = getdetJdV(cellvalues, q_point)
        x = spatial_coordinate(cellvalues, q_point, coords)
        ## Loop over test shape functions
        for i in 1:n_basefuncs
            δu  = shape_value(cellvalues, q_point, i)
            ∇δu = shape_gradient(cellvalues, q_point, i)
            ## Add contribution to fe
            fe[i] += analytical_rhs(x) * δu * dΩ
            ## Loop over trial shape functions
            for j in 1:n_basefuncs
                ∇u = shape_gradient(cellvalues, q_point, j)
                ## Add contribution to Ke
                Ke[i, j] += (∇δu ⋅ ∇u) * dΩ
            end
        end
    end
    return Ke, fe
end

# Standard assembly copy pasta for Poisson problem
function assemble_global(cellvalues::CellScalarValues, K::SparseMatrixCSC, dh::DofHandler)
    ## Allocate the element stiffness matrix and element force vector
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    ## Allocate global force vector f
    f = zeros(ndofs(dh))
    ## Create an assembler
    assembler = start_assemble(K, f)
    ## Loop over all cels
    for cell in CellIterator(dh)
        ## Reinitialize cellvalues for this cell
        reinit!(cellvalues, cell)
        coords = getcoordinates(cell)
        ## Compute element contribution
        assemble_element!(Ke, fe, cellvalues, coords)
        ## Assemble Ke and fe into K and f
        assemble!(assembler, celldofs(cell), Ke, fe)
    end
    return K, f
end

# Assemble and solve
function solve(dh, ch, cellvalues)
    K, f = assemble_global(cellvalues, create_sparsity_pattern(dh), dh);
    apply!(K, f, ch)
    u = K \ f;
end

# Check L2 convergence
function check_and_compute_convergence(dh, u, cellvalues; testatol)
    L2norm = 0.0
    L∞norm = 0.0
    local_test_passed = true
    for cell in CellIterator(dh)
        reinit!(cellvalues, cell)
        n_basefuncs = getnbasefunctions(cellvalues)
        coords = getcoordinates(cell)
        uₑ = u[celldofs(cell)]
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            x = spatial_coordinate(cellvalues, q_point, coords)
            uₐₙₐ    = prod(cos, x*π/2)
            uₐₚₚᵣₒₓ = function_value(cellvalues, q_point, uₑ)
            L2norm += norm(uₐₙₐ-uₐₚₚᵣₒₓ)*dΩ
            L∞norm = max(L∞norm, norm(uₐₙₐ-uₐₚₚᵣₒₓ))
            if local_test_passed && !isapprox(uₐₙₐ, uₐₚₚᵣₒₓ; atol=testatol)
                local_test_passed = false
                # Potentially print something
            end
        end
    end
    L2norm, L∞norm, local_test_passed
end

function check_solution(case::TestCase)
    dh, ch, cellvalues = setup_problem(case)
    u = solve(dh, ch, cellvalues)
    L2norm, L∞norm, local_test_passed = check_and_compute_convergence(dh, u, cellvalues; testatol=case.point_tol)
    L2_passed = (L2norm < case.L2_tol)
    L∞_passed = (L∞norm < case.L∞_tol)
    L2_passed || @show L2norm
    L∞_passed || @show L∞norm
    # Potentially print warnings...
    return (local_test_passed, L2_passed, L∞_passed)
end


# Finite element approximation parameters
cases = Dict(
    "Wedge, linear ip" => TestCase(element=Wedge, ip_order=1, qe_order=3, num_elems=25),
    "Triangle, linear ip" => TestCase(element=Triangle, ip_order=1, qe_order=2, num_elems=25))

for (name, case) in cases
    results = check_solution(case)
    if all(results)
        printstyled("$name: Passed \n"; color=:green, bold=true)
    else
        printstyled("$name: Failed: $results \n", color=:red, bold=true)
    end
end

[termi-official]

Looking at the current CI times for the examples in the documentation I do not think github runners have enough power to make such integration tests. We will likely just bump into a timeout before finishing such tests.

I had nothing specific in mind. With such scripts we can catch major regressions and actually compute convergence curves towards analytical solutions, which is super useful to see if everything works correctly together (because only then the convergence curve is correct).

[KnutAM]

But for the check, I suppose a single point (or maybe a few) on the curve would suffice (i.e. as I do in the modified script) - it would catch any regression in accuracy. And then this should be quite fast (I think docs etc. is taking a lot extra due to additional packages, but this case would only need Ferrite and its dependencies). Then for manual testing, one can still use the same functions to create the full convergence curves?
Or do you believe we need to check for quite fine meshes as well?

[termi-official]

I am not convinced that a single point is better than the regression tests which we already have in the examples. However, it can be deceiving, because faulty interpolations can still give acceptable results until we go deep down into the asymptotic convergence regime.

[termi-official]

This new test increases CI times by ~30 seconds. We can drive down the timings by removing some of the interpolations from the tests, but I am not sure which.

[fredrikekre]

Mistake?

[termi-official]

Yes. Not sure how this happened.

[termi-official]

I use this regularly, so I decided to export it here. If you think that this requires more discussion then I will revert it and adjust the test.

[termi-official]

If we want a larger infrastructure for convergence checking we might want a dedicated CI queue.