Skip to content
Snippets Groups Projects
Select Git revision
  • 5f557378360da84eabe369cbe41c4533ec07c510
  • master default protected
  • parallelistation_of_subdomain_loop
  • parallelistation_test
  • nonzero_gli_term_for_nonwetting_in_TPR
  • testing_updating_pwi_as_pwiminus1_into_calculation_of_pnwi
  • v2.2
  • v2.1
  • v2.0
  • v1.0
10 results

TP-TP-layered_soil_pure_dd.py

Blame
  • problems.jl 1.96 KiB
    using LinearAlgebra: UniformScaling, Diagonal, opnorm
    using StaticKernels: ExtensionNothing
    
    #"min_u 1/2 * |Au-f|_2^2 + λ*TV(u) + β/2 |u|_2 + γ |u|_1"
    "min_p 1/2 * |p₂ - div(p₁) - g|_B^2 + χ_{|p₁|≤λ} + χ_{|p₂|≤γ}"
    struct DualTVL1ROFOpProblem{Tg,TB,Tλ,Tγ} <: Problem
        "dual data"
        g::Tg
        "B operator"
        B::TB
        "TV regularization parameter"
        λ::Tλ
        "L1 regularization parameter"
        γ::Tγ
    end
    
    DualTVL1ROFOpProblem(g, B, λ) = DualTVL1ROFOpProblem(g, B, λ, nothing)
    DualTVL1ROFOpProblem(g, B, λ::Real) = DualTVL1ROFOpProblem(g, B, fill!(similar(g, typeof(λ)), λ))
    
    function energy(p, prob::DualTVL1ROFOpProblem)
        d = ndims(p)
    
        @inline kfΛ(w) = @inbounds divergence(w) + prob.g[w.position]
        kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
    
        # v = div(p) + g
        v = map(kΛ, extend(p, ExtensionNothing()))
    
        # |v|_B^2
        u = prob.B * vec(v)
        return sum(dot.(u, vec(v)))
    end
    
    function recover_u(p, prob::DualTVL1ROFOpProblem)
        d = ndims(prob.g)
        u = similar(prob.g)
        v = similar(prob.g)
    
        @inline kfΛ(w) = @inbounds divergence(w)
        kΛ = Kernel{ntuple(_->-1:1, d)}(kfΛ)
    
        # v = div(p) + A'*f
        map!(kΛ, v, StaticKernels.extend(p, StaticKernels.ExtensionNothing())) # extension: nothing
        v .+= prob.g
        # u = B * v
        mul!(vec(u), prob.B, vec(v))
        return u
    end
    
    mynorm(v) = norm(vec(v))
    @inline mynorm(v::SArray{<:Any,<:SArray{<:Any,Float64,1},1}) =
        sqrt(sum(sum(v[i] .^ 2) for i in eachindex(v)))
    
    function residual(p, prob::DualTVL1ROFOpProblem)
        d = ndims(p)
        grad = Kernel{ntuple(_->0:1, d)}(gradient)
    
        u = recover_u(p, prob)
        q = map(grad, StaticKernels.extend(u, StaticKernels.ExtensionReplicate()))
    
        res = q .- mynorm.(q) .* p ./ prob.λ
        return sum(dot.(res, res)) / length(p)
    end
    
    # operator norm of B
    normB(p::DualTVL1ROFOpProblem{<:Any,<:UniformScaling}) = p.B.λ
    normB(p::DualTVL1ROFOpProblem{<:Any,<:Diagonal}) = maximum(opnorm.(diag(p.B)))
    normB(p::DualTVL1ROFOpProblem) = opnorm(p.B)