430: Parameter Derivatives (stationary)

(source code)

Explore different ways to calculate sensitivities. This is still experimental.

module Example430_ParameterDerivativesStationary

using VoronoiFVM, ExtendableGrids
using GridVisualize
using ExtendableSparse
using ForwardDiff, DiffResults
using SparseDiffTools, SparseArrays
using ILUZero, LinearSolve

"""
    f(P)

Parameter dependent function which creates system and solves it
"""
function f(P; n = 10)
    p = P[1]

    valuetype = typeof(p)
    nspecies = 1
    ispec = 1

    function flux!(f, u, edge)
        f[1] = (1 + p) * (u[1, 1]^2 - u[1, 2]^2)
    end

    function r!(f, u, edge)
        f[1] = p * u[1]^5
    end

    function bc!(f, u, node)
        boundary_dirichlet!(f, u, node, ispec, 1, 0.0)
        boundary_dirichlet!(f, u, node, ispec, 3, p)
    end

    X = collect(0:(1.0 / n):1)
    grid = simplexgrid(X, X)
    sys = VoronoiFVM.System(grid; valuetype, species = [1], flux = flux!, reaction = r!,
                            bcondition = bc!)
    tff = VoronoiFVM.TestFunctionFactory(sys)
    tfc = testfunction(tff, [1], [3])
    sol = solve(sys; inival = 0.5)
    [integrate(sys, tfc, sol)[1]]
end

"""
    runf(;Plotter, n=10)

Run parameter series, plot f(p), df(p).
For each p,create a new system. Use VoronoiFVM with dual numbers. Pass parameters via closure.
"""
function runf(; Plotter = nothing, n = 10)
    P = 0.1:0.05:2
    dresult = DiffResults.JacobianResult(ones(1))
    F = zeros(0)
    DF = zeros(0)
    ff(p) = f(p; n)
    @time for p ∈ P
        ForwardDiff.jacobian!(dresult, ff, [p])
        push!(F, DiffResults.value(dresult)[1])
        push!(DF, DiffResults.jacobian(dresult)[1])
    end
    vis = GridVisualizer(; Plotter, legend = :lt)
    scalarplot!(vis, P, F; color = :red, label = "f")
    scalarplot!(vis, P, DF; color = :blue, label = "df", clear = false, show = true)
    sum(DF)
end

function fluxg!(f, u, edge, data)
    f[1] = (1 + data.p) * (u[1, 1]^2 - u[1, 2]^2)
end

function rg!(f, u, edge, data)
    f[1] = data.p * u[1]^5
end

function bcg!(f, u, node, data)
    boundary_dirichlet!(f, u, node, 1, 1, 0.0)
    boundary_dirichlet!(f, u, node, 1, 3, data.p)
end

Base.@kwdef mutable struct MyData{Tv}
    p::Tv = 1.0
end

"""
    rung(;Plotter, n=10)

Same as runf, but keep one system pass parameters via data.
"""
function rung(; Plotter = nothing, method_linear = SparspakFactorization(), n = 10)
    X = collect(0:(1.0 / n):1)
    grid = simplexgrid(X, X)

ugly but simple. By KISS we should first provide this way.

    sys = nothing
    data = nothing
    tfc = nothing

    function g(P)
        Tv = eltype(P)
        if isnothing(sys)
            data = MyData(one(Tv))
            sys = VoronoiFVM.System(grid; valuetype = Tv, species = [1], flux = fluxg!,
                                    reaction = rg!, bcondition = bcg!, data,
                                    unknown_storage = :dense)
            tff = VoronoiFVM.TestFunctionFactory(sys)
            tfc = testfunction(tff, [1], [3])
        end
        data.p = P[1]
        sol = solve(sys; inival = 0.5, method_linear, precon_linear = ILUZeroPreconditioner())
        [integrate(sys, tfc, sol)[1]]
    end

    dresult = DiffResults.JacobianResult(ones(1))

    P = 0.1:0.05:2
    G = zeros(0)
    DG = zeros(0)
    @time for p ∈ P
        ForwardDiff.jacobian!(dresult, g, [p])
        push!(G, DiffResults.value(dresult)[1])
        push!(DG, DiffResults.jacobian(dresult)[1])
    end

    vis = GridVisualizer(; Plotter, legend = :lt)
    scalarplot!(vis, P, G; color = :red, label = "g")
    scalarplot!(vis, P, DG; color = :blue, label = "dg", clear = false, show = true)
    sum(DG)
end

#########################################################################

function fluxh!(f, u, edge)
    p = parameters(u)[1]
    f[1] = (1 + p) * (u[1, 1]^2 - u[1, 2]^2)
end

function rh!(f, u, edge)
    p = parameters(u)[1]
    f[1] = p * u[1]^5
end

function bch!(f, u, node)
    p = parameters(u)[1]
    boundary_dirichlet!(f, u, node, 1, 1, 0.0)
    boundary_dirichlet!(f, u, node, 1, 3, p)
end

"""
    runh(;Plotter, n=10)

Same as runf, but use "normal" calculation (don't solve in dual numbers), and calculate dudp during
main assembly loop.

This needs quite a bit of additional implementation + corresponding API and still lacks local assembly of the
measurement derivative (when using testfunction based calculation) when calculating current.
"""
function runh(; Plotter = nothing, n = 10)
    X = collect(0:(1.0 / n):1)
    grid = simplexgrid(X, X)

    sys = VoronoiFVM.System(grid; species = [1], flux = fluxh!, reaction = rh!,
                            bcondition = bch!, unknown_storage = :dense, nparams = 1)
    tff = VoronoiFVM.TestFunctionFactory(sys)
    tfc = testfunction(tff, [1], [3])

    function measp(params, u)
        Tp = eltype(params)
        up = Tp.(u)
        integrate(sys, tfc, up; params = params)[1]
    end

    params = [0.0]

    function mymeas!(meas, U)
        u = reshape(U, sys)
        meas[1] = integrate(sys, tfc, u; params)[1]
        nothing
    end

    dp = 0.05
    P = 0.1:dp:2
    U0 = solve(sys; inival = 0.5, params = [P[1]])

    ndof = num_dof(sys)
    colptr = [i for i = 1:(ndof + 1)]
    rowval = [1 for i = 1:ndof]
    nzval = [1.0 for in = 1:ndof]
    ∂m∂u = zeros(1, ndof)
    colors = matrix_colors(∂m∂u)

    H = zeros(0)
    DH = zeros(0)
    DHx = zeros(0)
    m = zeros(1)

    @time for p ∈ P
        params[1] = p
        sol = solve(sys; inival = 0.5, params)
        mymeas!(m, sol)
        push!(H, m[1])

this one is expensive - we would need to assemble this jacobian via local calls

        forwarddiff_color_jacobian!(∂m∂u, mymeas!, vec(sol); colorvec = colors)

need to have the full derivative of m vs p

        ∂m∂p = ForwardDiff.gradient(p -> measp(p, sol), params)

        dudp = sys.matrix \ vec(sys.dudp[1])
        dmdp = -∂m∂u * dudp + ∂m∂p
        push!(DH, dmdp[1])
    end

    vis = GridVisualizer(; Plotter, legend = :lt)
    scalarplot!(vis, P, H; color = :red, label = "h")
    scalarplot!(vis, P, DH; color = :blue, label = "dh", clear = false, show = true)
    sum(DH)
end

using Test
function runtests()
    testval = 489.3432830184927
    @test runf()  ≈ testval
    @test rung()  ≈ testval
    @test runh()  ≈ testval
end

end

This page was generated using Literate.jl.