Fix DETPFlash reduced Gibbs scaling for GammaPhi/PTFlashWrapper (#517)

* Fix DETPFlash reduced Gibbs scaling for GammaPhi/PTFlashWrapper

- Fix DETPFlash objective scaling: accumulate phase reduced Gibbs (G/RT) and normalize by total moles to return molar reduced Gibbs g = G/(nRT) in Obj_de_tp_flash.
  - Previously: GammaPhi/PTFlashWrapper path produced reduced Gibbs, but the objective applied an extra /RT, causing inconsistent g vs RRTPFlash/MichelsenTPFlash.
- Make __eval_G_DETPFlash(model::EoSModel, ...) return reduced Gibbs consistently (g/RT) to match PTFlashWrapper/GammaPhi path.
- Fix gammaphi_gibbs empty-phase return: always return (g, v) tuple for sum(w)==0 to avoid type/unwrap inconsistency.
- Refactor gammaphi_gibbs typing: introduce a single promoted eltype and reuse it for both g and v zeros.
- Add minimal SASS ask–tell usage example at top of sass_at.jl.

	modified:   src/methods/property_solvers/multicomponent/tp_flash/DifferentialEvolutiontp_flash.jl
	modified:   src/models/CompositeModel/GammaPhi.jl
	modified:   src/modules/solvers/metaheuristics/sass_at.jl
	modified:   test/test_methods_api_flash.jl

* fix CI

Author: abcdvvvv

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

@@ -169,31 +169,33 @@ function Obj_de_tp_flash(model,p,T,n,dividers,numphases,x,nvals,vcache,logspace
     #nvals = zeros(TT,numphases, numspecies)
     #Calculate partition of species into phases
     partition!(dividers,n,x,nvals)
-    #Calculate Overall Gibbs energy (J)
+    #Accumulate total reduced Gibbs energy, G/(RT)
     #If any errors are encountered, return a big number, ensuring point is discarded
     #by DE Algorithm
     G = _0
     for i ∈ 1:numphases
         ni = @view(nvals[i, :])
-        gi,vi = __eval_G_DETPFlash(model,p,T,ni,equilibrium)
+        gi_reduced,vi = __eval_G_DETPFlash(model,p,T,ni,equilibrium)
         vcache[i] = vi
-        G += gi
+        G += gi_reduced
         #calling with PTn calls the internal volume solver
         #if it returns an error, is a bug in our part.
     end
     if logspace
         dividers .= log.(dividers)
     end
-    R̄ = Rgas(model)
-    return ifelse(isnan(G),bignum,G/R̄/T)
+    # Per the FlashData definition, return the molar reduced Gibbs energy (g = G/(nRT))
+    ∑n = sum(n)
+    return ifelse(isnan(G),bignum,G/∑n)
 end
 
 #indirection to allow overloading this evaluation in activity models
 function __eval_G_DETPFlash(model::EoSModel,p,T,ni,equilibrium)
     phase = is_lle(equilibrium) ? :liquid : :unknown
+    RT = Rgas(model)*T
     vi = volume(model,p,T,ni;phase = phase)
     g = VT_gibbs_free_energy(model, vi, T, ni)
-    return g,vi
+    return g/RT,vi
 end
 
 numphases(method::DETPFlash) = method.numphases

src/models/CompositeModel/GammaPhi.jl

@@ -279,11 +279,12 @@ end
 
 function gammaphi_gibbs(wrapper::PTFlashWrapper,p,T,w,phase = :unknown,vol = NaN)
     model = wrapper.model
+    TT = Base.promote_eltype(__γ_unwrap(model),p,T,w)
     RT = Rgas(model)*T
     ∑w = sum(w)
-    iszero(∑w) && return zero(Base.promote_eltype(model,p,T,w))
-    g_ideal = sum(xlogx,w) - xlogx(sum(w)) 
-    vl = zero(Base.promote_eltype(__γ_unwrap(model),p,T,w))
+    iszero(∑w) && return zero(TT), zero(TT)
+    g_ideal = sum(xlogx,w) - xlogx(∑w)
+    vl = zero(TT)
     if is_liquid(phase)
         return excess_gibbs_free_energy(__γ_unwrap(model),p,T,w)/RT + g_ideal,vl
     elseif is_vapour(phase)

src/modules/solvers/metaheuristics/sass_at.jl

@@ -1,3 +1,26 @@
+"""
+Simple usage example:
+
+```julia
+using Clapeyron
+const SOL = Clapeyron.Solvers
+
+f(x) = sum(abs2, x)
+dim = 2
+lb = fill(-5.0, dim)
+ub = fill(5.0, dim)
+
+algo = SOL.SASS(50, 10_000, lb, ub; seed = 1, stagnation_evals = 0, stagnation_tol = 0.0)
+while !SOL.isdone(algo)
+    x = SOL.ask!(algo)
+    y = f(x)
+    SOL.tell!(algo, y)
+end
+
+best_x, best_y = SOL.best(algo)
+```
+"""
+
 mutable struct Arch
     max_npop::Int
     x::Matrix{Float64}

test/test_methods_api_flash.jl

@@ -35,7 +35,7 @@
         #https://julialang.zulipchat.com/#narrow/channel/265161-Clapeyron.2Ejl/topic/The.20meaning.20of.20subcooled.20liquid.20flash.20results/near/534216551
         model_zulip2 = PR(["n-butane", "n-pentane", "n-hexane", "n-heptane"])
         res3 = Clapeyron.tp_flash2(model_zulip2, 1e5 , 450, z_zulip1, RRTPFlash(equilibrium=:vle))
-        
+
         if Clapeyron.numphases(res3) == 2
             @test isone(res3.fractions[2])
             @test res3.volumes[1] ≈ 0.03683358805181434 rtol = 1e-6
@@ -69,15 +69,15 @@
 
     @testset "DE Algorithm" begin
         #VLLE eq
-        @test Clapeyron.tp_flash(system, p, T, z, DETPFlash(numphases = 3))[3] ≈ -6.759674475174073 rtol = 1e-6
+        @test Clapeyron.tp_flash(system, p, T, z, DETPFlash(numphases = 3))[3] ≈ -6.759674475174073 rtol = 1e-12
         #LLE eq with activities
         act_system = UNIFAC(["water","cyclohexane","propane"])
         flash0 = Clapeyron.tp_flash(act_system, p, T, [0.5,0.5,0.0], DETPFlash(equilibrium = :lle))
         act_x0 = activity_coefficient(act_system, p, T, flash0[1][1,:]) .* flash0[1][1,:]
         act_y0 = activity_coefficient(act_system, p, T, flash0[1][2,:]) .* flash0[1][2,:]
-        @test Clapeyron.dnorm(act_x0,act_y0) < 1e-6
+        @test Clapeyron.dnorm(act_x0,act_y0) < 1e-5
         flash_RR = Clapeyron.tp_flash(act_system, p, T, [0.5,0.5,0.0], RRTPFlash(equilibrium = :lle))
-        @test flash0[3] ≈ flash_RR[3] rtol = 1e-12
+        @test flash0[3] ≈ flash_RR[3] rtol = 1e-11
     end
 
     @testset "Multiphase algorithm" begin
@@ -245,7 +245,7 @@
         flash4 = tp_flash(model_vle, 2500.0, 300.15, [0.9, 0.1], MichelsenTPFlash())
 
         @test flash4[1] ≈
-        [0.9239684120579815 0.07603158794201849; 
+        [0.9239684120579815 0.07603158794201849;
         0.793479931206839 0.20652006879316098] rtol = 1e-6
         #test equality of activities does not make sense in VLE
     end
@@ -488,7 +488,7 @@ end
     h475 = Clapeyron.PS.enthalpy(fluid475,1.742722525216547e6,-89.04935789018991,[1.0,1.0])
     res475 = Clapeyron.PS.flash(fluid475,1.742722525216547e6,-89.04935789018991,[1.0,1.0])
     @test enthalpy(fluid475,res475) ≈ h475 rtol = 1e-6
-    
+
     #issue 492
     fluid492 = GERG2008(["propane", "butane"])
     p492 = 1.5e6
@@ -503,7 +503,7 @@ end
     p_in_506  = 101325.0; z_506 = [1.0,1.0]
     T_in_506  = dew_temperature(fluid506,p_in_506,z_506)[1] + 10;  # we are now pure vapour
     p_out_506 = 3.0*p_in_506
-    h_in_506  = enthalpy(fluid506,p_in_506,T_in_506,z) 
+    h_in_506  = enthalpy(fluid506,p_in_506,T_in_506,z)
     s_in_506  = Clapeyron.PH.entropy(fluid506, p_in_506, h_in_506,z_506,phase = phase506)
     h_out_506 = Clapeyron.PS.enthalpy(fluid506, p_out_506, s_in_506,z_506,phase = phase506)
     s_out_506 = Clapeyron.PH.entropy(fluid506, p_out_506, h_out_506,z_506,phase = phase506)