add pereira composition generator to tpd

Author: longemen3000

src/methods/property_solvers/stability/tpd.jl

@@ -327,7 +327,7 @@ function _tpd(model,p,T,z,cache = tpd_cache(model,p,T,z),break_first = false,lle
 
     #create a list of what test composition strategies we use
     tpd_strategies = tpd_plan(model,z,isliquidz,lle,id_test,K_test,pure_test)
-
+    verbose && @info "bulk composition d(z) vector: $dz"
     for strategy in tpd_strategies
         w,phasew,skip = tpd_test_composition!(strategy,cond,w_test,K,dz,verbose)
         skip && continue
@@ -382,6 +382,8 @@ function tpd_plan(model,z,is_liquidz,lle,id_test,K_test,pure_test)
             end
         end
     end
+
+
     return plan
 end
 
@@ -404,8 +406,8 @@ function tpd_test_composition!(strategy,conds,w_test,K,dz,verbose)
         w_test ./= sum(w_test)
     elseif plan == :pure
         z_pure!(w_test,ix)
-    elseif plan == :electrolyte_balanced
-        z_electrolyte_balanced!(model,w_test,z,ixx)
+    elseif plan == :pereira
+        z_pereira!(w_test,z,ixx)
     elseif is_k_plan && !skip_k
         if all(iszero,K)
             K .= tp_flash_K0(model,p,T,z)
@@ -455,8 +457,8 @@ function tpd_print_strategy(strategy)
 
     elseif plan == :pure
         res = "Strategy: pure initial point, test phase: $phase"
-    elseif plan == :electrolyte_balanced
-        res = "Strategy: electrolyte-constrainted composition"
+    elseif plan == :pereira
+        res = "Strategy: Pereira composition generator"
     else
         res = ""
     end
@@ -504,19 +506,36 @@ function z_pure!(K,i)
     K
 end
 
-function z_electrolyte_balanced!(model,w_test,z,ixx)
-    isolv,ielec,ic = ixx
+function z_pereira!(w,z,ixx)
+    i,b,_ = ixx
+	n = length(z)
+    @assert n != i
 
-    Z = model.charge
-    Zx = Z[ic]
-    w_test .= 0
-    for i in 1:length(model)
-        isolv == i && (w_test[i] = z[i])
-        ielec == i && (w_test[i] = z[i])
+
+	lb = eps(eltype(w))*1e2
+	ub = 1.0 - lb
+
+    for j = 1:n
+        if (j == i)
+            if b > 0
+                w[j] = 0.5*z[i]
+            else
+                w[j] = z[i] + 0.5*(1.0 - z[i])
+            end
+        else
+            if b > 0
+                w[j] = (1.0 - 0.5*z[i])/(n-1)
+            else
+                w[j] = (1.0 - (z[i] + 0.5*(1.0 - z[i])))/(n-1)
+            end
+        end
     end
-    Zw = dot(Z,w_test)
-    w_test[ic] = - Zw/Zx
-    w_test ./= sum(w_test)
+
+    w .= clamp.(w,lb,ub)
+    sumw = sum(w)
+    w ./= sumw
+    return w
+
 end
 
 function z_norm(z,w)
@@ -637,7 +656,6 @@ function K0_lle_init(model::EoSModel, p, T, z)
     return K
 end
 
-
 """
 
 

src/models/Electrolytes/ISElectrolyte.jl

@@ -51,10 +51,11 @@ end
 Rgas(model::ISElectrolyteWrapper) = Rgas(model.model)
 Rgas(model::ISElectrolyteIdealWrapper) = Rgas(model.model)
 
+#=
 function eos_impl(model::ISElectrolyteWrapper,V,T,z)
     w = to_ion(model.salt,z)
     return eos_impl(model.model,V,T,w)
-end
+end=#
 
 function tp_flash_K0!(K,model::ISElectrolyteModel,p,T,z)
     neutral = ones(Bool,length(model))
@@ -73,8 +74,6 @@ function each_split_model(model::ISElectrolyteWrapper,I_salt)
     return ISElectrolyteWrapper(model.components[I_salt],each_split_model(model.model,I_ion),salt_i)
 end
 
-export ISElectrolyteWrapper
-
 function volume_impl(model::ISElectrolyteWrapper, p, T, z, phase, threaded, vol0)
     w = to_ion(model.salt,z)
     return volume(model.model, p, T, w, phase=phase, threaded=threaded, vol0=vol0)
@@ -97,4 +96,6 @@ end
 function T_scale(model::ISElectrolyteWrapper,z)
     w = to_ion(model.salt,z)
     return T_scale(model.model,w)
-end
+end
+
+export ISElectrolyteWrapper

src/models/Electrolytes/stability.jl

@@ -1,4 +1,49 @@
-function tpd_plan(model::ElectrolyteModel,z,is_liquidz,lle,id_test,K_test,pure_test)
+function tpd_plan(model::ISElectrolyteModel,z,is_liquidz,lle,id_test,K_test,pure_test)
+    plan = Tuple{Symbol,Symbol,NTuple{3,Int}}[]
+    nc = length(model)
+    neutral = ones(Bool,length(model))
+    isalts = model.salt.isalts
+    neutral[isalts] .= false
+
+    if is_liquidz && id_test && !lle && iszero(length(isalts))
+        push!(plan,(:ideal_gas,:vapour,0))
+    end
+
+    ids = sortperm(z)
+    if is_liquidz
+        for i in 1:nc
+            #idx_solvent = ids[i]
+            #if neutral[idx_solvent]
+            #    push!(plan,(:pure,:liquid,(idx_solvent,0,0)))
+            #end
+            if i != nc
+                push!(plan,(:pereira,:liquid,(i,-1,0)))
+                push!(plan,(:pereira,:liquid,(i,1,0)))
+            end
+        end
+        if !lle
+            for i in 1:nc
+                !iszero(ZZ[i]) && push!(plan,(:pure,:vapour,(ids[i],0,0)))
+            end
+        end
+    else
+        for i in 1:nc
+            #idx_solvent = ids[i]
+            #if neutral[idx_solvent]          
+            #    push!(plan,(:pure,:liquid,(idx_solvent,0,0)))
+            #end
+
+            if i != nc
+                push!(plan,(:pereira,:liquid,(i,-1,0)))
+                push!(plan,(:pereira,:liquid,(i,1,0)))
+            end
+        end
+    end
+    return plan
+end
+
+#=
+function tpd_plan(model::ESElectrolyteModel,z,is_liquidz,lle,id_test,K_test,pure_test)
     plan = Tuple{Symbol,Symbol,NTuple{3,Int}}[]
 
     if is_liquidz && id_test && !lle
@@ -59,4 +104,5 @@ function tpd_plan(model::ElectrolyteModel,z,is_liquidz,lle,id_test,K_test,pure_t
         end
     end
     return plan
-end
+end
+=#