module yaeos__equilibria_boundaries_phase_envelopes_mp !! Multiphase PT envelope calculation module. !! !! This module contains the functions to calculate the PT envelope of a !! mixture with multiple phases. use yaeos__constants, only: pr, R use yaeos__equilibria_equilibrium_state, only: EquilibriumState use yaeos__equilibria_boundaries_auxiliar, only: & detect_critical, check_critical_jump use yaeos__models_ar, only: ArModel use yaeos__math, only: solve_system implicit none private public :: PTEnvelMP public :: pt_F_NP public :: pt_envelope type :: PTEnvelMP !! Multiphase PT envelope. type(MPPoint), allocatable :: points(:) !! Array of converged points. real(pr), allocatable :: Tc(:) !! Critical temperatures. real(pr), allocatable :: Pc(:) !! Critical pressures. contains procedure :: write => write_envelope_PT_MP procedure, nopass :: solve_point procedure, nopass :: get_values_from_X end type PTEnvelMP type :: MPPoint !! Multiphase equilibria point. integer :: np !! Number of phases integer :: nc !! Number of components real(pr) :: beta_w !! Fraction of the reference (incipient) phase. real(pr), allocatable :: betas(:) !! Fractions of the main phases. real(pr) :: P !! Pressure [bar] real(pr) :: T !! Temperature [K] real(pr), allocatable :: x_l(:, :) !! Mole fractions of the main phases. real(pr), allocatable :: w(:) !! Mole fractions of the incipient phase. character(len=14), allocatable :: kinds_x(:) !! Kinds of the main phases. character(len=14) :: kind_w !! Kind of the reference phase. integer :: iters !! Number of iterations needed to converge the point. integer :: ns !! Number of the specified variable. end type MPPoint contains type(PTEnvelMP) function pt_envelope(& model, z, np, kinds_x, kind_w, x_l0, w0, betas0, P0, T0, ns0, dS0, beta_w, points, & max_pressure & ) !! # `pt_envelope` !! Calculation of a multiphase PT envelope. !! !! # Description !! Calculates a PT envelope is calculated using the continuation method. !! The envelope is calculated by solving the system of equations for each !! point of the envelope. The system of equations is solved using the !! Newton-Raphson method. !! !! This function requires the system specification conditions, which are !! the fluid composition (\z\), the number of phases that are not !! incipient; defined as \(np\), proper intialization values, the !! variables that end with `0` are the initial guess; the mole fraction !! of the reference phase `beta_w` which when it is equal to 0 means that !! we are calculating a phase boundary. use yaeos__auxiliar, only: optval class(ArModel), intent(in) :: model real(pr), intent(in) :: z(:) !! Mixture global composition. integer, intent(in) :: np !! Number of main phases. character(len=14), intent(in) :: kinds_x(np) !! Kind of the main phases. character(len=14), intent(in) :: kind_w !! Kind of the reference phase. real(pr), intent(in) :: x_l0(np, size(z)) !! Initial guess for the mole fractions of each phase. arranged as !! an array of size `(np, nc)`, where nc is the number of components !! and `np` the number of main phases. Each row correspond to the !! composition of each main phaase. real(pr), intent(in) :: w0(size(z)) !! Initial guess for the mole fractions of the !! reference/incipient phase. real(pr), intent(in) :: betas0(np) !! Initial guess for the fractions of the main phases. arranged as !! an array of size `(np)`, where `np` is the number of main phases. real(pr), intent(in) :: P0 !! Initial guess for the pressure [bar]. real(pr), intent(in) :: T0 !! Initial guess for the temperature [K]. integer, intent(in) :: ns0 !! Number of the specified variable. !! The variable to be specified. This is the variable that will be !! used to calculate the first point of the envelope. The variable !! can be any of the variables in the vector X, but it is recommended !! to use the temperature or pressure. The variables are aranged as !! follows: !! !! - `X(1:nc*np) = ln(K_i^l)`: \(\frac{x_i^l}{w_i}\) !! - `X(nc*np+1:nc*np+np) = \beta_i^l`: Fraction of each main phase. !! - `X(nc*np+np+1) = ln(P)`: Pressure [bar]. !! - `X(nc*np+np+2) = ln(T)`: Temperature [K]. real(pr), intent(in) :: dS0 !! Step size of the specification for the next point. !! This is the step size that will be used to calculate the next point. !! Inside the algorithm this value is modified to adapt the step size !! to facilitate the convergence of each point. real(pr), intent(in) :: beta_w !! Fraction of the reference (incipient) phase. integer, optional, intent(in) :: points !! Number of points to calculate. real(pr), optional, intent(in) :: max_pressure !! Maximum pressure [bar] to calculate. !! If the pressure of the point is greater than this value, the !! calculation is stopped. !! This is useful to avoid calculating envelopes that go to infinite !! values of pressure. type(MPPoint), allocatable :: env_points(:) !! Array of converged points. type(MPPoint) :: point !! Converged point. real(pr) :: max_P !! Maximum pressure [bar] to calculate. real(pr) :: F(size(z) * np + np + 2) !! Vector of functions valuated. real(pr) :: dF(size(z) * np + np + 2, size(z) * np + np + 2) !! Jacobian matrix. real(pr) :: dXdS(size(z) * np + np + 2) !! Sensitivity of the variables wrt the specification. real(pr) :: X(size(z) * np + np + 2) !! Vector of variables. real(pr) :: dX(size(z) * np + np + 2) !! Step for next point estimation. integer :: nc !! Number of components. integer :: its !! Number of iterations to solve the current point. integer :: max_iterations = 100 !! Maximum number of iterations to solve the point. integer :: number_of_points !! Number of points to calculate. real(pr) :: x_l(np, size(z)) !! Mole fractions of the main phases. real(pr) :: w(size(z)) !! Mole fractions of the incipient phase. real(pr) :: betas(np) !! Fractions of the main phases. real(pr) :: P !! Pressure [bar]. real(pr) :: T !! Temperature [K]. integer :: i !! Point calculation index integer :: iT !! Index of the temperature variable. integer :: iP !! Index of the pressure variable. integer :: lb !! Lower bound, index of the first component of a phase integer :: ub !! Upper bound, index of the last component of a phase integer :: inner !! Number of times a failed point is retried to converge integer :: ns !! Number of the specified variable real(pr) :: dS !! Step size of the specification for the next point real(pr) :: S !! Specified value real(pr) :: X0(size(X)) !! Initial guess for the point real(pr) :: X_last_converged(size(X)) !! Last converged point real(pr) :: Xc(size(X)) !! Vector of variables at the critical point logical :: found_critical !! If true, a critical point was found logical :: jumped_critical !! If true, a critical point was jumped character(len=14) :: x_kinds(np), w_kind real(pr) :: Tc !! Critical temperature [K] real(pr) :: Pc !! Critical pressure [bar] real(pr) :: Vl(np) !! Molar volumes of the main phases [L/mol] real(pr) :: Vw !! Molar volume of the reference phase [L/mol] integer :: l !! Phase index nc = size(z) iP = np*nc + np + 1 iT = np*nc + np + 2 number_of_points = optval(points, 1000) max_P = optval(max_pressure, 2000._pr) do i=1,np lb = (i-1)*nc + 1 ub = i*nc X(lb:ub) = log(x_l0(i, :)/w0) end do X(np*nc + 1:np*nc + np) = betas0 X(np*nc + np + 1) = log(P0) X(np*nc + np + 2) = log(T0) ns = ns0 S = X(ns) dS = dS0 x_kinds = kinds_x w_kind = kind_w allocate(env_points(0), pt_envelope%Tc(0), pt_envelope%Pc(0)) F = 1 its = 0 X0 = X call solve_point(& model, z, np, beta_w, x_kinds, w_kind, X, ns, S, dXdS, & F, dF, Vl, Vw, its, 1000 & ) do i=1,number_of_points X0 = X call solve_point(& model, z, np, beta_w, x_kinds, w_kind, X, ns, S, dXdS, & F, dF, Vl, Vw, its, max_iterations & ) ! The point might not converge, in this case we try again with an ! initial guess closer to the previous (converged) point. inner = 0 do while(i > 1 .and. its >= max_iterations .and. inner < 10) inner = inner + 1 X = X0 - (1 - real(inner, pr) / 10._pr) * dX if (ns > 0) then S = X(ns) end if call solve_point(& model, z, np, beta_w, x_kinds, w_kind, X, ns, S, dXdS, & F, dF, Vl, Vw, its, max_iterations & ) end do ! Convert the values of the vector of variables into human-friendly ! variables. call get_values_from_X(X, np, z, beta_w, x_l, w, betas, P, T) ! Attach the new point to the envelope. point = MPPoint(& np=np, nc=nc, betas=betas, P=P, T=T, x_l=x_l, w=w, beta_w=beta_w, & kinds_x=x_kinds, kind_w=w_kind, iters=its, ns=ns & ) ! Check if the system is close to a critical point, and try to jump ! over it. call check_critical_jump(nc, np, ns, X, X_last_converged, Xc, jumped_critical) if (jumped_critical) then Tc = exp(Xc(iT)) Pc = exp(Xc(iP)) pt_envelope%Tc = [pt_envelope%Tc, Tc] pt_envelope%Pc = [pt_envelope%Pc, Pc] end if ! detect_critical jumps over the CP, so we save the X_last_converged ! before that X_last_converged = X call detect_critical(& nc, np, i, x_kinds, w_kind, .false., & X_last_converged, X, dXdS, ns, dS, S, & found_critical, Xc) ! Update the specification for the next point. call update_specification(its, nc, np, X, dF, dXdS, ns, S, dS, Vl, Vw) ! If the point did not converge, stop the calculation if (& any(isnan(F)) .or. its > max_iterations & .or. P > max_P & .or. any(betas < -1e-14) .or. any(betas > 1 + 1e-14) & .or. abs(dS) <= 1e-14 & .or. (T < 100._pr .and. P < 1e-2_pr) & ) exit env_points = [env_points, point] ! Next point estimation. dX = dXdS * dS do while(abs(exp(X(iT)) - exp(X(iT) + dX(iT))) > 7) dX = dX/2 end do do while(abs(exp(X(iP)) - exp(X(iP) + dX(iP))) > 5) dX = dX/2 end do X = X + dX if (ns > 0) then S = X(ns) end if end do ! This moves the locally saved points to the output variable. call move_alloc(env_points, pt_envelope%points) end function pt_envelope subroutine pt_F_NP(model, z, np, beta_w, kinds_x, kind_w, X, ns, S, F, dF, Vl, Vw) !! Function to solve at each point of a multi-phase envelope. use iso_fortran_env, only: error_unit class(ArModel), intent(in) :: model !! Model to use. real(pr), intent(in) :: z(:) !! Mixture global composition. integer, intent(in) :: np !! Number of main phases. real(pr), intent(in) :: beta_w !! Fraction of the reference (incipient) phase. character(len=14), intent(in) :: kinds_x(np) !! Kind of the main phases. character(len=14), intent(in) :: kind_w !! Kind of the reference phase. real(pr), intent(in) :: X(:) !! Vector of variables. integer, intent(in) :: ns !! Number of specification. real(pr), intent(in) :: S !! Specification value. real(pr), intent(out) :: F(size(X)) !! Vector of functions valuated. real(pr), intent(out) :: df(size(X), size(X)) !! Jacobian matrix. real(pr), intent(out) :: Vw real(pr), intent(out) :: Vl(np) ! X variables real(pr) :: K(np, size(z)) real(pr) :: P real(pr) :: T real(pr) :: betas(np) ! Main phases variables real(pr) :: moles(size(z)) real(pr), dimension(np, size(z)) :: x_l, lnphi_l, dlnphi_dt_l, dlnphi_dp_l real(pr), dimension(np, size(z), size(z)) :: dlnphi_dn_l real(pr) :: dpdv_l(np), dpdn_l(np, size(z)), dpdt_l(np), dPdT real(pr) :: lnphi(size(z)), dlnphi_dt(size(z)), dlnphi_dp(size(z)) real(pr), dimension(size(z), size(z)) :: dlnphi_dn real(pr) :: dpdv_w, dpdn_w(size(z)), dpdT_w real(pr) :: dPdV, dPdN(size(z)) real(pr) :: dVwdn(size(z)), dVwdT, dVwdP ! Incipient phase variables real(pr), dimension(size(z)) :: w, lnphi_w, dlnphi_dt_w, dlnphi_dp_w real(pr), dimension(size(z), size(z)) :: dlnphi_dn_w ! Derivatives of w wrt beta and K real(pr) :: dwdb(np, size(z)) real(pr) :: dwdlnK(np, size(z)) real(pr) :: denom(size(z)) real(pr) :: denomdlnK(np, size(z), size(z)) real(pr) :: dx_l_dlnK(np, np, size(z)) integer :: i, j, l, phase, nc integer :: lb, ub integer :: idx_1, idx_2 integer :: nv !! Specified volume for ln(V(nv)/Vw) nc = size(z) ! ======================================================================== ! Extract variables from the vector X ! ------------------------------------------------------------------------ do l=1,np lb = (l-1)*nc + 1 ub = l*nc K(l, :) = exp(X(lb:ub)) end do betas = X(np*nc + 1:np*nc + np) P = exp(X(np*nc + np + 1)) T = exp(X(np*nc + np + 2)) denom = 0 denom = matmul(betas, K) + beta_w denomdlnK = 0 do i=1,nc denomdlnK(:, i, i) = betas(:)*K(:, i) end do w = z/denom ! ======================================================================== ! Calculation of fugacities coeficients and their derivatives ! ------------------------------------------------------------------------ call model%lnphi_pt(& w, P, T, V=Vw, root_type=kind_w, lnphi=lnphi_w, & dlnphidp=dlnphi_dp_w, dlnphidt=dlnphi_dt_w, dlnphidn=dlnphi_dn_w, & dPdV=dpdv_w, dPdN=dpdn_w, dPdT=dPdT_w & ) do l=1,np x_l(l, :) = K(l, :)*w call model%lnphi_pt(& x_l(l, :), P, T, V=Vl(l), root_type=kinds_x(l), lnphi=lnphi, & dlnphidp=dlnphi_dp, dlnphidt=dlnphi_dt, dlnphidn=dlnphi_dn, & dPdV=dpdv, dPdN=dPdN, dPdT=dpdt & ) lnphi_l(l, :) = lnphi dlnphi_dn_l(l, :, :) = dlnphi_dn dlnphi_dt_l(l, :) = dlnphi_dt dlnphi_dp_l(l, :) = dlnphi_dp dpdv_l(l) = dpdv dpdn_l(l, :) = dPdN dpdt_l(l) = dpdt end do ! ======================================================================== ! Calculation of the system of equations ! ------------------------------------------------------------------------ do l=1,np ! Select the limits of the function lb = (l-1)*nc + 1 ub = l*nc F(lb:ub) = X(lb:ub) + lnphi_l(l, :) - lnphi_w F(nc * np + l) = sum(x_l(l, :) - w) end do F(nc * np + np + 1) = sum(betas) + beta_w - 1 if (ns < 0) then nv = abs(ns) F(nc * np + np + 2) = log(Vl(nv)/Vw) - S else F(nc * np + np + 2) = X(ns) - S endif ! ======================================================================== ! Derivatives and Jacobian Matrix of the whole system ! ------------------------------------------------------------------------ df = 0 dwdlnK = 0 do l=1,np ! Save the derivatives of w wrt beta and K of the incipient phase dwdb(l, :) = -z * K(l, :)/denom**2 dwdlnK(l, :) = -K(l, :) * betas(l)*z/denom**2 end do do l=1,np do phase=1,np dx_l_dlnK(phase, l, :) = dwdlnK(l, :) * K(phase, :) if (phase == l) then dx_l_dlnK(phase, l, :) = dx_l_dlnK(phase, l, :) + w * K(l, :) end if end do end do do l=1,np ! Derivatives of the isofugacity equations ! wrt lnK do phase=1,np do i=1, nc do j=1,nc idx_1 = i + (phase-1)*nc idx_2 = j + (l-1)*nc df(idx_1, idx_2) = & dlnphi_dn_l(phase, i, j) * dx_l_dlnK(phase, l, j) & - dlnphi_dn_w(i, j) * dwdlnK(l, j) if (i == j .and. phase == l) then df(idx_1, idx_2) = df(idx_1, idx_2) + 1 end if end do end do end do ! wrt betas do j=1,np lb = (j-1)*nc + 1 ub = j*nc do i=1,nc df(lb+i-1, np*nc + l) = & sum(K(j, :) * dlnphi_dn_l(j, i, :)*dwdb(l, :) & - dlnphi_dn_w(i, :)*dwdb(l, :)) end do end do ! wrt T,p do i=1,nc lb = (l-1)*nc + i df(lb, nc*np+np+1) = P*(dlnphi_dp_l(l, i) - dlnphi_dp_w(i)) df(lb, nc*np+np+2) = T*(dlnphi_dt_l(l, i) - dlnphi_dt_w(i)) end do ! ==================================================================== ! Derivatives of the sum of mole fractions ! -------------------------------------------------------------------- ! wrt lnK do phase=1,np do j=1,nc lb = nc*np + phase ub = j + (l-1)*nc df(lb, ub) = df(lb, ub) + (dx_l_dlnK(phase, l, j) - dwdlnK(l, j)) end do end do ! wrt beta do j=1,np lb = nc*np + j df(lb,np*nc+l) = sum(K(j, :) * dwdb(l, :) - dwdb(l, :)) end do ! Derivatives of sum(beta)==1 df(nc * np + np + 1, np*nc + l) = 1 end do block real(pr) :: dVwdlnKl(np, size(z)) real(pr) :: dVnvdlnKl(np, size(z)) real(pr) :: dVwdb(np) real(pr) :: dVnvdb(np) real(pr) :: dVnvdn(size(z)), dVnvdP, dVnvdT real(pr) :: Vnv integer :: spec_eq if (ns < 0) then spec_eq = nc * np + np + 2 nv = abs(ns) Vnv = Vl(nv) dVwdn = dpdn_w / dpdv_w dVwdP = 1.0_pr / dpdv_w dVwdT = dPdT_w / dpdv_w dVnvdn = dpdn_l(nv, :) / dpdv_l(nv) dVnvdP = 1.0_pr / dpdv_l(nv) dVnvdT = dPdT_l(nv) / dpdv_l(nv) do l=1,np dVwdlnKl(l, :) = dVwdn * dwdlnK(l, :) dVnvdlnKl(l, :) = dVnvdn * dx_l_dlnK(nv, l, :) dVwdb(l) = sum(dVwdn * dwdb(l, :)) dVnvdb(l) = sum(dVnvdn * K(nv, :) * dwdb(l, :)) end do do l=1,np lb = (l-1) * nc ! wrt lnK^l_i do i=1,nc df(spec_eq, lb + i) = & dVwdlnKl(l, i) / Vw - dVnvdlnKl(l, i) / Vnv end do ! wrt beta lb = np * nc + l df(spec_eq, lb) = & dVwdb(l) / Vw - dVnvdb(l) / Vnv end do df(spec_eq, nc * np + np + 1) = Vw / Vnv * P * (1 / Vw * dVnvdP - Vnv / Vw**2 * dVwdP) df(spec_eq, nc * np + np + 2) = -Vw / Vnv * T * (1 / Vw * dVnvdT - Vnv / Vw**2 * dVwdT) else df(nc * np + np + 2, ns) = 1 end if end block end subroutine pt_F_NP subroutine solve_point(model, z, np, beta_w, kinds_x, kind_w, X, ns, S, dXdS, F, dF, Vl, Vw, iters, max_iterations) use iso_fortran_env, only: error_unit use yaeos__math, only: solve_system use stdlib_linalg, only: eye class(ArModel), intent(in) :: model !! Model to use. real(pr), intent(in) :: z(:) !! Mixture global composition. integer, intent(in) :: np !! Number of main phases real(pr), intent(in) :: beta_w !! Fraction of the reference (incipient) phase character(len=14), intent(in) :: kinds_x(np) !! Kind of the main phases character(len=14), intent(in) :: kind_w !! Kind of the reference phase real(pr), intent(in out) :: X(:) !! Vector of variables integer, intent(in) :: ns !! Number of specification real(pr), intent(in) :: S !! Specification value real(pr), intent(in) :: dXdS(size(X)) real(pr), intent(out) :: F(size(X)) !! Vector of functions valuated real(pr), intent(out) :: df(size(X), size(X)) !! Jacobian matrix real(pr), intent(out) :: Vw !! Reference phase volume real(pr), intent(out) :: Vl(np) !! Main phases volumes integer, intent(in) :: max_iterations !! Maximum number of iterations to solve the point integer, intent(out) :: iters !! Number of iterations to solve the current point integer :: i integer :: l !! Phase index integer :: iT !! Temperature index integer :: iP !! Pressure index integer :: iBetas(np) !! Indices of the betas in X integer :: nc !! Number of components real(pr) :: X0(size(X)) real(pr) :: dX(size(X)) ! Homotopy solving real(pr) :: H(size(X)), dH(size(X), size(X)) real(pr) :: G(size(X)), dG(size(X), size(X)) real(pr) :: t logical :: can_solve nc = size(z) iP = np*nc + np + 1 iT = np*nc + np + 2 X0 = X can_solve = .true. iBetas = [(i, i=np*nc+1, np*nc+np)] do iters=1,max_iterations call pt_F_NP(model, z, np, beta_w, kinds_x, kind_w, X, ns, S, F, dF, Vl, Vw) if (any(isnan(F)) .and. can_solve) then X = X - 0.9 * dX can_solve = .false. cycle end if dX = solve_system(dF, -F) if (maxval(abs(F)) < 1e-9_pr) exit X = X + dX if (iters == 10) then ! If after 10 iterations it is not converging try fixed-point ! homotopy. Homotopy consist of solving a modified version of the ! system of equations that is easier to solve, and gradually ! transforming it into the original system. ! The function to solve is defined as: ! ! \[ ! H = t F + (1-t) G ! \] ! Where \( G = X - X_0 \) is a simple function that has a known ! solution at \( X = X_0 \), and \( t \) ! goes from 0 to 1 in 10 steps. X0 = X H = 10 do i=0,10 do while(maxval(abs(H)) > 1e-9_pr) t = real(i, pr)/10._pr call pt_F_NP(model, z, np, beta_w, kinds_x, kind_w, X, ns, S, F, dF, Vl, Vw) G = (X - X0) dG = eye(size(X)) H = t * F + (1-t) * G dH = t * dF + (1-t) * dG dX = solve_system(dH, -H) X = X + dX end do end do end if end do end subroutine solve_point subroutine update_specification(its, nc, np, X, dF, dXdS, ns, S, dS, Vl, Vw) !! # update_specification !! Change the specified variable for the next step. !! !! # Description !! Using the information of a converged point and the Jacobian matrix of !! the function. It is possible to determine the sensitivity of the !! variables with respect to the specification. This information is used !! to update the specification for the next point. Choosing the variable !! with the highest sensitivity. !! This can be done by solving the system of equations: !! !! \[ !! J \frac{dX}{dS} + \frac{dF}{dS} = 0 !! \] !! !! for the \( \frac{dX}{dS} \) vector. The variable with the highest value !! of \( \frac{dX}{dS} \) is chosen as the new specification. !! !! # References !! integer, intent(in) :: its !! Iterations to solve the current point. integer, intent(in) :: nc !! Number of components in the mixture. integer, intent(in) :: np !! Number of main phases. real(pr), intent(in out) :: X(:) !! Vector of variables. real(pr), intent(in out) :: dF(:, :) !! Jacobian matrix. real(pr), intent(in out) :: dXdS(:) !! Sensitivity of the variables wrt the specification. integer, intent(in out) :: ns !! Number of the specified variable. real(pr), intent(in out) :: S !! Specified value. real(pr), intent(in out) :: dS !! Step size of the specification for the next point. real(pr), intent(in) :: Vl(:) !! Molar volumes of the main phases [L/mol]. real(pr), intent(in) :: Vw !! Molar volume of the reference phase [L/mol]. real(pr) :: dFdS(size(X)) !! Sensitivity of the functions wrt the specification. integer :: i integer :: l !! Phase index integer :: lb !! Lower bound of each phase integer :: ub !! Upper bound of each phase integer :: iT integer :: iP integer :: iBetas(np) real(pr) :: dT, dP iBetas = [(i, i=np*nc+1, np*nc+np)] iP = size(X) - 1 iT = size(X) dFdS = 0 dFdS(size(X)) = -1 dXdS = solve_system(dF, -dFdS) ns = maxloc(abs(dXdS), dim=1) ! ======================================================================== ! For each phase, check if the mole fractions are too low. ! this can be related to criticality and it is useful to force the ! specification of compositions. ! ------------------------------------------------------------------------ dS = dXdS(ns)*dS dXdS = dXdS/dXdS(ns) dS = dS * 3./its dS = sign(min(abs(dS), 0.1_pr), dS) do l=1,np lb = (l-1)*nc + 1 ub = l*nc if (maxval(abs(X(lb:ub))) < 0.2) then ns = maxloc(abs(dXdS(lb:ub)), dim=1) + lb - 1 dS = dXdS(ns)*dS dXdS = dXdS/dXdS(ns) dS = sign(0.01_pr, dS) exit end if end do end subroutine update_specification subroutine get_values_from_X(X, np, z, beta_w, x_l, w, betas, P, T) !! # get_values_from_X !! Extract the values of the variables from the vector X. !! real(pr), intent(in) :: X(:) !! Vector of variables. integer, intent(in) :: np !! Number of main phases. real(pr), intent(in) :: z(:) !! Mixture composition. real(pr), intent(in) :: beta_w !! Fraction of the reference phase. real(pr), intent(out) :: x_l(np, size(z)) !! Mole fractions of the main phases. real(pr), intent(out) :: w(size(z)) !! Mole fractions of the incipient phase. real(pr), intent(out) :: betas(np) !! Fractions of the main phases. real(pr), intent(out) :: P !! Pressure [bar]. real(pr), intent(out) :: T !! Temperature [K]. integer :: nc !! Number of components. integer :: i !! Loop index. integer :: lb !! Lower bound of each phase. integer :: ub !! Upper bound of each phase. nc = size(z) betas = X(np*nc + 1:np*nc + np) P = exp(X(np*nc + np + 1)) T = exp(X(np*nc + np + 2)) ! Extract the K values from the vector of variables do i=1,np lb = (i-1)*nc + 1 ub = i*nc x_l(i, :) = exp(X(lb:ub)) end do ! Calculate the mole fractions of the incipient phase w = z/(matmul(betas, x_l) + beta_w) ! Calculate the mole fractions of the main phases with the previously ! calculated K values do i=1,np x_l(i, :) = x_l(i, :)*w end do end subroutine get_values_from_X subroutine write_envelope_PT_MP(env, unit) class(PTEnvelMP), intent(in) :: env integer, intent(in) :: unit integer :: i, j integer :: np, nc real(pr) :: P, T real(pr), allocatable :: betas(:) real(pr), allocatable :: w(:) real(pr), allocatable :: x_l(:, :) np = size(env%points) nc = size(env%points(1)%w) do i=1,np P = env%points(i)%P T = env%points(i)%T betas = env%points(i)%betas w = env%points(i)%w x_l = env%points(i)%x_l write(unit, "(*(E15.5,2x))") P, T, betas, w, (x_l(j, :), j=1, size(x_l,dim=1)) end do end subroutine write_envelope_PT_MP end module yaeos__equilibria_boundaries_phase_envelopes_mp