430: Parameter Derivatives (stationary)
Explore different ways to calculate sensitivities. This is still experimental.
module Example430_ParameterDerivativesStationaryusing VoronoiFVM, ExtendableGridsusing GridVisualizeusing ExtendableSparseusing ExtendableSparse: ILUZeroPreconBuilderusing ForwardDiff, DiffResultsusing SparseArraysusing ILUZero, LinearSolveusing DifferentiationInterface, SparseConnectivityTracer, SparseMatrixColorings""" f(P)Parameter dependent function which creates system and solves it"""function f(P; n = 10, method_linear = SparspakFactorization()) p = P[1] valuetype = typeof(p) nspecies = 1 ispec = 1 function flux!(f, u, edge, data) f[1] = (1 + p) * (u[1, 1]^2 - u[1, 2]^2) return nothing end function r!(f, u, edge, data) f[1] = p * u[1]^5 return nothing end function bc!(f, u, node, data) boundary_dirichlet!(f, u, node, ispec, 1, 0.0) boundary_dirichlet!(f, u, node, ispec, 3, p) return nothing 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, method_linear) return [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, method_linear = SparspakFactorization()) P = 0.1:0.05:2 dresult = DiffResults.JacobianResult(ones(1)) F = zeros(0) DF = zeros(0) ff(p) = f(p; n, method_linear) @time for p in 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) return sum(DF)endfunction fluxg!(f, u, edge, data) f[1] = (1 + data.p) * (u[1, 1]^2 - u[1, 2]^2) return nothingendfunction rg!(f, u, edge, data) f[1] = data.p * u[1]^5 return nothingendfunction bcg!(f, u, node, data) boundary_dirichlet!(f, u, node, 1, 1, 0.0) boundary_dirichlet!(f, u, node, 1, 3, data.p) return nothingendBase.@kwdef mutable struct MyData{Tv} p::Tv = 1.0end""" rung(;Plotter, n=10)Same as runf, but keep one system pass parameters via data."""function rung(; Plotter = nothing, n = 10, method_linear = SparspakFactorization()) 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) return [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 in 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) return sum(DG)end#########################################################################function fluxh!(f, u, edge, data) p = parameters(u)[1] f[1] = (1 + p) * (u[1, 1]^2 - u[1, 2]^2) return nothingendfunction rh!(f, u, edge, data) p = parameters(u)[1] f[1] = p * u[1]^5 return nothingendfunction bch!(f, u, node, data) p = parameters(u)[1] boundary_dirichlet!(f, u, node, 1, 1, 0.0) boundary_dirichlet!(f, u, node, 1, 3, p) return nothingend""" runh(;Plotter, n=10)Same as runf, but use "normal" calculation (don't solve in dual numbers), and calculate dudp duringmain assembly loop.This needs quite a bit of additional implementation + corresponding API and still lacks local assembly of themeasurement derivative (when using testfunction based calculation) when calculating current."""function runh(; Plotter = nothing, n = 10, method_linear = SparspakFactorization()) 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) return 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] return nothing end dp = 0.05 P = 0.1:dp:2 state = VoronoiFVM.SystemState(sys) U0 = solve!(state; inival = 0.5, params = [P[1]]) input = VoronoiFVM.dofs(U0) output = zeros(1) backend = AutoSparse( AutoForwardDiff(); sparsity_detector = TracerSparsityDetector(), coloring_algorithm = GreedyColoringAlgorithm() ) jac_prep = prepare_jacobian(mymeas!, output, backend, input; strict = Val(false)) ∂m∂u = similar(sparsity_pattern(jac_prep), Float64) H = zeros(0) DH = zeros(0) DHx = zeros(0) m = zeros(1) @time for p in P params[1] = p sol = solve!(state; inival = 0.5, params, method_linear) mymeas!(m, sol) push!(H, m[1])this one is expensive - we would need to assemble this jacobian via local calls
DifferentiationInterface.jacobian!(mymeas!, output, ∂m∂u, jac_prep, backend, vec(sol))need to have the full derivative of m vs p
∂m∂p = ForwardDiff.gradient(p -> measp(p, sol), params) dudp = state.matrix \ vec(state.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) return sum(DH)endusing Testfunction runtests() testval = 489.3432830184927 for method_linear in [UMFPACKFactorization(), SparspakFactorization()] @test runf(; method_linear) ≈ testval @test rung(; method_linear) ≈ testval @test runh(; method_linear) ≈ testval end return nothingendendThis page was generated using Literate.jl.