001: New linear solver API

(source code)

module Example001_Solvers

# under development

using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize
using LinearSolve
using ExtendableSparse
using AMGCLWrap
using AlgebraicMultigrid
using LinearAlgebra
using Test

function main(; n = 10, Plotter = nothing, assembly = :edgwwise, kwargs...)
    h = 1.0 / convert(Float64, n)
    X = collect(0.0:h:1.0)
    Y = collect(0.0:h:1.0)

    grid = simplexgrid(X, Y)
    nn = num_nodes(grid)

    eps = 1.0e-2

    function reaction(f, u, node, data)
        f[1] = u[1]^2
        return nothing
    end

    function flux(f, u, edge, data)
        f[1] = eps * (u[1, 1]^2 - u[1, 2]^2)
        return nothing
    end

    function source(f, node, data)
        x1 = node[1] - 0.5
        x2 = node[2] - 0.5
        f[1] = exp(-20.0 * (x1^2 + x2^2))
        return nothing
    end

    function storage(f, u, node, data)
        f[1] = u[1]
        return nothing
    end

    function bcondition(f, u, node, data)
        boundary_dirichlet!(
            f,
            u,
            node;
            species = 1,
            region = 2,
            value = ramp(node.time; dt = (0, 0.1), du = (0, 1))
        )
        boundary_dirichlet!(
            f,
            u,
            node;
            species = 1,
            region = 4,
            value = ramp(node.time; dt = (0, 0.1), du = (0, 1))
        )
        return nothing
    end

    sys = VoronoiFVM.System(
        grid; reaction, flux, source, storage, bcondition, assembly,
        species = [1]
    )
    @info "UMFPACK:"
    umf_sol = solve(sys; inival = 0.5, method_linear = UMFPACKFactorization(), kwargs...)

    @info "KLU:"
    sol = solve(sys; inival = 0.5, method_linear = KLUFactorization(), kwargs...)
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Sparspak:"
    sol = solve(sys; inival = 0.5, method_linear = SparspakFactorization(), kwargs...)
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov-ilu0:"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = ILUZeroPreconditioner(),
        kwargs...
    )
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov-block1"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = BlockPreconditioner(;
            partitioning = [1:(nn ÷ 2), (nn ÷ 2 + 1):nn],
            factorization = ILU0Preconditioner()
        ),
        kwargs...
    )
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov-block2"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = BlockPreconditioner(;
            partitioning = [1:2:nn, 2:2:nn],
            factorization = UMFPACKFactorization()
        ),
        kwargs...
    )
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov - delayed factorization:"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = SparspakFactorization(),
        keepcurrent_linear = false,
        kwargs...
    )
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov - jacobi:"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = JacobiPreconditioner(),
        keepcurrent_linear = true,
        kwargs...
    )
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov - SA_AMG:"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = SA_AMGPreconditioner(),
        keepcurrent_linear = true,
        kwargs...
    )
    @test norm(sol - umf_sol, Inf) < 1.0e-7

    @info "Krylov - AMGCL_AMG:"
    sol = solve(
        sys;
        inival = 0.5,
        method_linear = KrylovJL_BICGSTAB(),
        precon_linear = AMGCL_AMGPreconditioner(),
        keepcurrent_linear = true,
        kwargs...
    )
    return @test norm(sol - umf_sol, Inf) < 1.0e-7

end

function runtests()
    @testset "edgewise" begin
        main(; assembly = :edgewise)
    end
    @testset "cellwise" begin
        main(; assembly = :cellwise)
    end
    return nothing
end
end

This page was generated using Literate.jl.