NestedADError in MultiPhaseTPFlash

[spinnau]

In Clapeyron 0.6.21 using MultiPhaseTPFlash results in NestedADError in some cases:

 using Clapeyron

components = ["carbon dioxide", "water"]
model = EOS_CG(components)

T = 293.15
p = 1e6
x = [0.5, 0.5]

res = Clapeyron.tp_flash2(model, p, T, x, MultiPhaseTPFlash())

results in the error

NestedADError: Found nested Duals of type ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f#multi_g_obj##0"{MultiFluid{EmpiricAncillary, AsymmetricMixing, EmpiricDeparture}, Float64, Float64, Vector{Float64}, Int64, Int64, Base.RefValue{Int64}}, Float64}, ForwardDiff.Dual{ForwardDiff.Tag{Clapeyron.var"#f#multi_g_obj##0"{MultiFluid{EmpiricAncillary, AsymmetricMixing, EmpiricDeparture}, Float64, Float64, Vector{Float64}, Int64, Int64, Base.RefValue{Int64}}, Float64}, Float64, 5}, 5}. This is currently not supported in implicit differentiation.

Stacktrace:
  [1] nested_ad_check(a::Vector{ForwardDiff.Dual{ForwardDiff.Tag{…}, ForwardDiff.Dual{…}, 5}})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl:283](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl#line=282)
  [2] (::Base.var"#foreach##0#foreach##1"{typeof(Clapeyron.nested_ad_check)})(::Nothing, x::Vector{ForwardDiff.Dual{…}})
    @ Base ./tuple.jl:691
  [3] BottomRF
    @ ./reduce.jl:78 [inlined]
  [4] afoldl(::Base.BottomRF{…}, ::Nothing, ::MultiFluid{…}, ::Float64, ::Float64, ::Vector{…})
    @ Base ./operators.jl:602
  [5] _foldl_impl(op::Base.BottomRF{…}, init::Nothing, itr::Tuple{…})
    @ Base ./reduce.jl:60
  [6] foldl_impl(op::Base.BottomRF{…}, nt::Nothing, itr::Tuple{…})
    @ Base ./reduce.jl:40
  [7] mapfoldl_impl(f::typeof(identity), op::Base.var"#foreach##0#foreach##1"{…}, nt::Nothing, itr::Tuple{…})
    @ Base ./reduce.jl:36
  [8] mapfoldl(f::Function, op::Function, itr::Tuple{MultiFluid{…}, Float64, Float64, Vector{…}}; init::Nothing)
    @ Base ./reduce.jl:167
  [9] mapfoldl
    @ ./reduce.jl:167 [inlined]
 [10] #foldl#271
    @ ./reduce.jl:190 [inlined]
 [11] foreach(f::typeof(Clapeyron.nested_ad_check), itr::Tuple{MultiFluid{…}, Float64, Float64, Vector{…}})
    @ Base ./tuple.jl:691
 [12] nested_ad_check(a::Tuple{MultiFluid{…}, Float64, Float64, Vector{…}})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl:288](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl#line=287)
 [13] implicit_ad_check(a::Tuple{MultiFluid{…}, Float64, Float64, Vector{…}})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl:320](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl#line=319)
 [14] __gradients_for_root_finders(x::StaticArraysCore.SVector{…}, tups::Tuple{…}, tups_primal::Tuple{…}, f::Clapeyron.var"#f_vec#__gradients_for_root_finders##0"{…})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl:229](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl#line=228)
 [15] __gradients_for_root_finders(x::Float64, tups::Tuple{…}, primal_tups::Tuple{…}, f::Clapeyron.var"#f#volume_ad##0")
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl:247](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/differentials.jl#line=246)
 [16] volume_ad(v::Float64, tups::Tuple{…}, tups_primal::Tuple{…})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/volume.jl:306](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/volume.jl#line=305)
 [17] _volume(model::MultiFluid{…}, p::Float64, T::Float64, z::Vector{…}, phase::Symbol, threaded::Bool, vol0::Nothing)
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/volume.jl:295](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/volume.jl#line=294)
 [18] volume(model::MultiFluid{…}, p::Float64, T::Float64, z::Vector{…}; phase::Symbol, threaded::Bool, vol0::Nothing)
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/volume.jl:273](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/volume.jl#line=272)
 [19] (::Clapeyron.var"#f#multi_g_obj##0"{…})(𝕏::Vector{…})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl:1162](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl#line=1161)
 [20] vector_mode_dual_eval!
    @ [~/.julia/packages/ForwardDiff/egQMG/src/apiutils.jl:24](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/apiutils.jl#line=23) [inlined]
 [21] vector_mode_gradient!(result::DiffResults.MutableDiffResult{…}, f::Clapeyron.var"#f#multi_g_obj##0"{…}, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff [~/.julia/packages/ForwardDiff/egQMG/src/gradient.jl:105](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/gradient.jl#line=104)
 [22] gradient!
    @ [~/.julia/packages/ForwardDiff/egQMG/src/gradient.jl:39](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/gradient.jl#line=38) [inlined]
 [23] InnerGradientForHess
    @ [~/.julia/packages/ForwardDiff/egQMG/src/hessian.jl:52](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/hessian.jl#line=51) [inlined]
 [24] vector_mode_dual_eval!(f!::ForwardDiff.InnerGradientForHess{…}, cfg::ForwardDiff.JacobianConfig{…}, y::Vector{…}, x::Vector{…})
    @ ForwardDiff [~/.julia/packages/ForwardDiff/egQMG/src/apiutils.jl:31](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/apiutils.jl#line=30)
 [25] vector_mode_jacobian!(result::Matrix{…}, f!::ForwardDiff.InnerGradientForHess{…}, y::Vector{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff [~/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl:157](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl#line=156)
 [26] jacobian!
    @ [~/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl:82](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/jacobian.jl#line=81) [inlined]
 [27] hessian!
    @ [~/.julia/packages/ForwardDiff/egQMG/src/hessian.jl:69](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/hessian.jl#line=68) [inlined]
 [28] hessian!
    @ [~/.julia/packages/ForwardDiff/egQMG/src/hessian.jl:67](http://localhost:8888/home/spinnau/.julia/packages/ForwardDiff/egQMG/src/hessian.jl#line=66) [inlined]
 [29] (::Clapeyron.Solvers.var"#fgh#ADScalarObjective##2"{…})(df::Vector{…}, d2f::Matrix{…}, x::Vector{…})
    @ Clapeyron.Solvers [~/.julia/packages/Clapeyron/UsPZJ/src/modules/solvers/optimize.jl:25](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/modules/solvers/optimize.jl#line=24)
 [30] upto_hessian
    @ [~/.julia/packages/NLSolvers/gB9BS/src/objectives.jl:68](http://localhost:8888/home/spinnau/.julia/packages/NLSolvers/gB9BS/src/objectives.jl#line=67) [inlined]
 [31] upto_hessian
    @ [~/.julia/packages/NLSolvers/gB9BS/src/optimize/problem_types.jl:26](http://localhost:8888/home/spinnau/.julia/packages/NLSolvers/gB9BS/src/optimize/problem_types.jl#line=25) [inlined]
 [32] init_f∇fB
    @ [~/.julia/packages/NLSolvers/gB9BS/src/quasinewton/approximations/newton.jl:19](http://localhost:8888/home/spinnau/.julia/packages/NLSolvers/gB9BS/src/quasinewton/approximations/newton.jl#line=18) [inlined]
 [33] prepare_variables(prob::NLSolvers.OptimizationProblem{…}, approach::NLSolvers.LineSearch{…}, x0::Vector{…}, ∇fz::Vector{…}, B::Nothing)
    @ NLSolvers [~/.julia/packages/NLSolvers/gB9BS/src/optimize/problem_types.jl:310](http://localhost:8888/home/spinnau/.julia/packages/NLSolvers/gB9BS/src/optimize/problem_types.jl#line=309)
 [34] _solve(mstyle::NLSolvers.InPlace, problem::NLSolvers.OptimizationProblem{…}, s0::Tuple{…}, approach::NLSolvers.LineSearch{…}, options::NLSolvers.OptimizationOptions{…}, cache::NLSolvers.QNVars{…})
    @ NLSolvers [~/.julia/packages/NLSolvers/gB9BS/src/optimize/linesearch/quasinewton.jl:61](http://localhost:8888/home/spinnau/.julia/packages/NLSolvers/gB9BS/src/optimize/linesearch/quasinewton.jl#line=60)
 [35] solve(problem::NLSolvers.OptimizationProblem{…}, x0::Vector{…}, approach::NLSolvers.LineSearch{…}, options::NLSolvers.OptimizationOptions{…})
    @ NLSolvers [~/.julia/packages/NLSolvers/gB9BS/src/optimize/linesearch/quasinewton.jl:26](http://localhost:8888/home/spinnau/.julia/packages/NLSolvers/gB9BS/src/optimize/linesearch/quasinewton.jl#line=25)
 [36] optimize(f::Function, x0::Vector{…}, method::NLSolvers.LineSearch{…}, options::NLSolvers.OptimizationOptions{…}; bounds::Nothing)
    @ Clapeyron.Solvers [~/.julia/packages/Clapeyron/UsPZJ/src/modules/solvers/optimize.jl:129](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/modules/solvers/optimize.jl#line=128)
 [37] optimize
    @ [~/.julia/packages/Clapeyron/UsPZJ/src/modules/solvers/optimize.jl:125](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/modules/solvers/optimize.jl#line=124) [inlined]
 [38] tp_flash_multi_neq!(model::MultiFluid{…}, p::Float64, T::Float64, z::Vector{…}, _result::Tuple{…}, ss_cache::Tuple{…}, options::MultiPhaseTPFlash{…})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl:1124](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl#line=1123)
 [39] tp_flash_multi(model::MultiFluid{…}, p::Float64, T::Float64, nn::Vector{…}, options::MultiPhaseTPFlash{…})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl:281](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl#line=280)
 [40] tp_flash_impl
    @ [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl:184](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl#line=183) [inlined]
 [41] tp_flash2(model::MultiFluid{…}, p::Float64, T::Float64, n::Vector{…}, method::MultiPhaseTPFlash{…})
    @ Clapeyron [~/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/flash/PT.jl:87](http://localhost:8888/home/spinnau/.julia/packages/Clapeyron/UsPZJ/src/methods/property_solvers/multicomponent/flash/PT.jl#line=86)
 [42] top-level scope
    @ In[9]:10
 [43] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
Some type information was truncated. Use `show(err)` to see complete types.

For other parameters like this changed pressure it works fine:

T = 293.15
p = 1e5
x = [0.5, 0.5]

Clapeyron.tp_flash2(model, p, T, x, MultiPhaseTPFlash()) 

Flash result at T = 293.15, p = 100000.0 with 2 phases:
 (x = [0.000878293, 0.999122], β = 0.488411, v = 1.80608e-5)
 (x = [0.976509, 0.0234905], β = 0.511589, v = 0.0242403)

[longemen3000]

Hmm, I know what is happening, the bug was introduced in the last release when unifying the solver to allow activity models, hopefully the fix should be simple

[longemen3000]

i got the system to work with MultiPhaseTPFlash and MichelsenTPFlash and added tests for the specific erroring case. Can you try the master branch, or do you need a new release?

[spinnau]

Thanks for looking into this. Current master fixes this error for me. I've tested with plotting some phase diagrams with a lot more calls of the flash routines.