reduce boxing variables (using JuliaLang/julia#60478)

Author: longemen3000

src/estimation/estimation.jl

@@ -357,16 +357,16 @@ function objective_function(estimation::Estimation,guesses)
         inputs = data.inputs
         outputs = data.outputs
         weights = data.weights
-        if isempty(inputs)
-            prediction = property(model_r)
+        prediction = if isempty(inputs)
+            property(model_r)
         elseif length(inputs)==1
-            prediction = property.(Ref(model_r),inputs[1])
+            property.(Ref(model_r),inputs[1])
         elseif length(inputs)==2
-            prediction = property.(Ref(model_r),inputs[1],inputs[2])
+            property.(Ref(model_r),inputs[1],inputs[2])
         elseif length(inputs)==3
-            prediction = property.(Ref(model_r),inputs[1],inputs[2],inputs[3])
+            property.(Ref(model_r),inputs[1],inputs[2],inputs[3])
         else
-            prediction = property.(Ref(model_r),inputs...)
+            property.(Ref(model_r),inputs...)
         end
 
         if length(outputs)==1

src/methods/initial_guess.jl

@@ -1015,9 +1015,9 @@ function solve_2ph_taylor(v10,v20,a1,da1,d2a1,a2,da2,d2a2,p_scale = 1.0,μ_scale
     end
     x0 = SVector((log(v10),log(v20)))
     x = Solvers.nlsolve2(F0,x0,Solvers.Newton2Var())
-    v1,v2 = exp(x[1]), exp(x[2])
-    p1 = log(v1/v10)*(-v1*d2a1) - da1
-    return v1, v2, p1
+    v1_sol,v2_sol = exp(x[1]), exp(x[2])
+    p1_sol = log(v1_sol/v10)*(-v1_sol*d2a1) - da1
+    return v1_sol, v2_sol, p1_sol
 end
 
 function solve_2ph_taylor(model1::EoSModel,model2::EoSModel,T,v1,v2,p_scale = 1.0,μ_scale = 1.0)

src/methods/property_solvers/fugacity_coefficient.jl

@@ -147,10 +147,11 @@ function ∂lnϕ∂n∂P(model::EoSModel, p, T, z=SA[1.], cache = ∂lnϕ_cache(
     n = sum(z)
     Z = p*V/RT/n
 
+    ncomponents = length(z)
     F_res(model, V, T, z) = eos_res(model, V, T, z) / RT
     fun(aux) = F_res(model, aux[1], T, @view(aux[2:(ncomponents+1)]))
 
-    ncomponents = length(z)
+    
     aux[1] = V
     aux[2:end] = z
     result = ForwardDiff.hessian!(result, fun, aux, hconfig, Val{false}())

src/methods/property_solvers/multicomponent/flash/QP.jl

@@ -46,8 +46,8 @@ function qp_flash_x0(model,β,p,z,method::FlashMethod)
             #Tmin,Tmax = extrema(T0)
             #we approximate sat(T) ≈ exp(-dpdT*T*T(1/T - 1/T0)/p)*p
             K = similar(dpdT,typeof(Tmax))
-            x = z ./ sum(z)
-            ft(_T) = qp_f0_T!(K,x,dpdT,p,_T,β)
+            xx = z ./ sum(z)
+            ft(_T) = qp_f0_T!(K,xx,dpdT,p,_T,β)
             #we do a search over Tmin-Tmax domain, finding the minimum value of the objective function
             Tm = β*Tmax + (1 - β)*Tmin
             Tr1 = range(Tmin,Tm,5*length(model))

src/methods/property_solvers/multicomponent/flash/QT.jl

@@ -42,8 +42,8 @@ function qt_flash_x0(model,β,T,z,method::FlashMethod)
             p_bubble = @sum(ps[i]*z[i])/∑z
             p_dew = ∑z/@sum(z[i]/ps[i])
             pmin,pmax = p_dew,p_bubble
-            x = z ./ sum(z)
-            fp(p) = qt_f0_p!(K,x,p,ps,β)
+            xx = z ./ sum(z)
+            fp(p) = qt_f0_p!(K,xx,p,ps,β)
             pm = β*pmin + (1-β)*pmax
             pr1 = range(pmin,pm,5*length(model))
             pr2 = range(pm,pmax,5*length(model))

src/methods/property_solvers/multicomponent/flash/general_flash.jl

@@ -573,17 +573,15 @@ function xy_flash(model::EoSModel,spec::FlashSpecifications,z,comps0,β0,volumes
     s = copy(input)
     #config,μconfig = xy_flash_config(model,input)
     f!(output,input) = xy_flash_neq(output,model,zz,np,input,new_spec,nothing)
-    Θ(_f) = 0.5*dot(_f,_f)
     srtol = abs2(cbrt(rtol))
     function Θ(_f,_z)
         f!(_f,_z)
         _f[slacks] .= 0
-        Θ(_f)
+        0.5*dot(_f,_f)
     end
     config = ForwardDiff.JacobianConfig(f!,F,x)
     #ForwardDiff.jacobian!(J,f!,F,x,config,Val{false}())
-    f!(F,x)
-    Θx = Θ(F)
+    Θx = Θ(F,x)
     Fnorm = sqrt(2*Θx)
     spec_norm = norm(viewlast(F,2),Inf)
     converged = Fnorm < rtol

src/methods/property_solvers/multicomponent/solids/sle_solubility.jl

@@ -31,9 +31,9 @@ function sle_solubility(model::CompositeModel,p,T,z;solute=nothing,x0=nothing)
         Tm = model.solid.params.Tm.values[idx_sol_s][1]
 
         idx_sol_l = zeros(Bool,nf)
-        solute_l = mapping[idx_sol_s][1]
-        ν_l = [solute_l[1][i][2] for i in 1:length(solute_l[1])]
-        solute_l = [solute_l[1][i][1] for i in 1:length(solute_l[1])]
+        solute_l1 = mapping[idx_sol_s][1]
+        ν_l = [solute_l1[1][i][2] for i in 1:length(solute_l1[1])]
+        solute_l = [solute_l1[1][i][1] for i in 1:length(solute_l1[1])]
 
 
         for i in solute_l
@@ -152,9 +152,9 @@ function sle_solubility_T(model::CompositeModel,z,p=1e5;solute=nothing,x0=nothin
         Tm = model.solid.params.Tm.values[idx_sol_s][1]
 
         idx_sol_l = zeros(Bool,nf)
-        solute_l = mapping[idx_sol_s][1]
-        ν_l = [solute_l[1][i][2] for i in 1:length(solute_l[1])]
-        solute_l = [solute_l[1][i][1] for i in 1:length(solute_l[1])]
+        solute_l1 = mapping[idx_sol_s][1]
+        ν_l = [solute_l1[1][i][2] for i in 1:length(solute_l1[1])]
+        solute_l = [solute_l1[1][i][1] for i in 1:length(solute_l1[1])]
 
         for i in solute_l
             idx_sol_l[fluid_components .== i] .= true

src/methods/property_solvers/multicomponent/tp_flash/electrolyte_flash.jl

@@ -283,7 +283,7 @@ end
 
 function rachfordrice(K, z, Z; β0=nothing, ψ0=nothing, non_inx=FillArrays.Fill(false,length(z)), non_iny=FillArrays.Fill(false,length(z)),verbose = false)
     # Function to solve Rachdord-Rice mass balance
-    β,status,limits = rachfordrice_β0(K.*exp.(Z.*ψ0),z,β0,non_inx,non_iny)
+    status = rachfordrice_status(K.*exp.(Z.*ψ0),z,β0,non_inx,non_iny)
     if status == RREq
         function rachford_rice_donnan(x,K,z,Z)
             β = x[1]
@@ -304,9 +304,9 @@ function rachfordrice(K, z, Z; β0=nothing, ψ0=nothing, non_inx=FillArrays.Fill
         ff(F,x) = rachford_rice_donnan(x,K,z,Z)
         results = Solvers.nlsolve(ff,x0)
         sol = Clapeyron.Solvers.x_sol(results)
-        β = sol[1]
-        ψ = sol[2]
-        return SVector(Base.promote(β,ψ))
+        β_sol = sol[1]
+        ψ_sol = sol[2]
+        return SVector(Base.promote(β_sol,ψ_sol))
     else
         return SVector(Base.promote(β0,ψ0))
     end

src/methods/property_solvers/multicomponent/tp_flash/multiphase.jl

@@ -1057,63 +1057,27 @@ function split_phase_tpd(model,p,T,z,w,phase_z = :unknown,phase_w = :unknown,vz
     x2 = similar(w)
     g1 = eos(model,vw,T,w) + p*vw
     gz = eos(model,vz,T,z) + p*vz
+    
     if isone(β2)
         x2 .= z
     else
         @. x2 = (z -  β2 * w) / (1 - β2)
     end
+
     x3 = x2
-    phase = phase_w == phase_z ? phase_z : :unknown
-    if is_vapour(phase_w) && is_unknown(phase)
-        phase = :liquid
-    end
-    v2 = volume(model,p,T,x2,threaded = false,phase = phase)
-    v3 = volume(model,p,T,x3,threaded = false,phase = phase)
-    g2 = eos(model,v2,T,x2) + p*v2
-    g3 = eos(model,v3,T,x3) + p*v3
-    dg1 = g1 - gz
-    dg2 = β2*g1 + (1-β2)*g2 - gz
-    function f(βx)
-        x3 .= (z .-  βx .* w) ./ (1 .- βx)
-        v3 = volume(model,p,T,x3,threaded = false,phase = phase)
-        g3 = eos(model,v3,T,x3) + p*v3
-        return βx*g1 + (1-βx)*g3 - gz
-    end
-    ϕ = 0.6180339887498949
-    βi = ϕ*β1 + (1-ϕ)*β2
-    βi0 = one(βi)*10
-    _1 = one(βi)
-    dgi = βi*g1 + (1-βi)*g3 - gz
-    βi*g1 + (1-βi)*g3 - gz
-    for i in 1:20
-        x3 .= (z .-  βi .* w) ./ (1 .- βi)
-        v3 = volume(model,p,T,x3,threaded = false,phase = phase)
-        g3 = eos(model,v3,T,x3) + p*v3
-        isnan(g3) && break
-        dgi = βi*g1 + (1-βi)*g3 - gz
-        #quadratic interpolation
-        A = @SMatrix [β1*β1 β1 _1; β2*β2 β2 _1; βi*βi βi _1]
-        b = SVector(dg1,dg2,dgi)
-        c = A\b
-        βi_intrp = -0.5*c[2]/c[1]
-        if dgi < dg2 < dg1
-            dg1,β1 = dgi,βi
-        elseif dgi < dg1 < dg2
-            dg2,β2,v2 = dgi,βi,v3
-        elseif dgi < dg1
-            dg1,β1,v1 = dgi,βi,v3
-        elseif dgi < dg2
-            dg2,β2,v2 = dgi,βi,v3
-        else
-            break
-        end
-        βi0 = βi
-        βi_bisec = ϕ*β1 + (1-ϕ)*β2
-        βi = β1 <= βi_intrp <= β2 ? βi_intrp : βi_bisec
-        abs(βi0 - βi) < 1e-5 && break
+    phase = phase_w == phase_z ? phase_z : (is_vapour(phase_w) ? :liquid : :unknown)
+    vcache = Ref(g1)
+
+    function ff(_β)
+        x3 .= (z .-  _β .* w) ./ (1 .- _β)
+        _v3  = volume(model,p,T,x3,threaded = false,phase = phase)
+        g3 = eos(model,_v3,T,x3) + p*_v3
+        vcache[] = _v3
+        return _β*g1 + (1-_β)*g3 - gz
     end
-    #@assert βi*w + (1-βi)*x3 ≈ z
-    return (1-βi),x3,v3,βi,w,vw,dgi
+
+    βi = Solvers.optimize(ff,(β2,oneunit(β2)),Solvers.BoundOptim1Var(),NLSolvers.OptimizationOptions(x_abstol = 1e-5))    #@assert βi*w + (1-βi)*x3 ≈ z
+    return (1-βi),x3,vcache[],βi,w,vw,dgi
 end
 
 function split_phase_k(model,p,T,z,K = nothing,vz = volume(model,p,T,z),pures = split_pure_model(model))

src/methods/property_solvers/singlecomponent/widom.jl

@@ -104,7 +104,7 @@ function widom_temperature(model,p;T0 = nothing,v0 = nothing)
         _T0 = T0*_1
     end
 
-    function f_wp(x::AbstractVector)
+    function f_wp_vec(x::AbstractVector)
         v,T = exp(x[1]),x[2]
         Cₚ(∂v,∂T) = VT_isobaric_heat_capacity(model,∂v,∂T,SA[1.0])
         fp(∂v,∂T) = pressure(model,∂v,∂T,SA[1.0])
@@ -119,12 +119,12 @@ function widom_temperature(model,p;T0 = nothing,v0 = nothing)
     #we do a refinement of the initial value first. it improves convergence near the critical point
     function f_wp(T)
         v = volume(model,p,T,phase = :l,vol0 = _v0)
-        return f_wp(SVector(log(v),T))[1]
+        return f_wp_vec(SVector(log(v),T))[1]
     end
     prob = Roots.ZeroProblem(f_wp,_T0)
     _T00 = Roots.solve(prob,Roots.Order2(),atol = 1e-1)
     _v00 = volume(model,p,_T00,phase = :l, vol0 = _v0)
-    res = Solvers.nlsolve2(f_wp,SVector(log(_v00),_T00),Solvers.Newton2Var())
+    res = Solvers.nlsolve2(f_wp_vec,SVector(log(_v00),_T00),Solvers.Newton2Var())
     log_v_widom,T_widom = res
     return T_widom,exp(log_v_widom)
 end
@@ -233,7 +233,7 @@ function ciic_temperature(model,p;T0 = nothing,v0 = nothing)
     lb_v = lb_volume(model,_T0,SA[1.0])
     scale = p*p/_T0/lb_volume(model,_T0,SA[1.0])
     #P*P/(T*v)
-    function f_ciic(x::AbstractVector)
+    function f_ciic_vec(x::AbstractVector)
         v,T = exp(x[1]),x[2]
         Cₚ(∂v) = VT_isobaric_heat_capacity(model,∂v,T,SA[1.0])
         fp(∂v) = pressure(model,∂v,T,SA[1.0])
@@ -247,13 +247,13 @@ function ciic_temperature(model,p;T0 = nothing,v0 = nothing)
     #we do a refinement of the initial value first. it improves convergence near the critical point
     function f_ciic(T)
         v = volume(model,p,T,phase = :l,vol0 = _v0)
-        return f_ciic(SVector(log(v),T))[1]
+        return f_ciic_vec(SVector(log(v),T))[1]
     end
     prob = Roots.ZeroProblem(f_ciic,_T0)
     _T00 = Roots.solve(prob,Roots.Order2(),atol = 1e-1)
     _v00 = volume(model,p,_T00,phase = :l, vol0 = _v0)
 
-    res = Solvers.nlsolve2(f_ciic,SVector(log(_v00),_T00),Solvers.Newton2Var())
+    res = Solvers.nlsolve2(f_ciic_vec,SVector(log(_v00),_T00),Solvers.Newton2Var())
     v_ciic,T_ciic = exp(res[1]),res[2]
     return T_ciic,v_ciic
 end