From f1e929350bf7a5329744f62478aca4ee3d9bd65a Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 7 Feb 2024 11:40:52 -0300 Subject: [PATCH] starting to write --- app/main.f90 | 215 +++++++---- fpm.toml | 5 +- src/Flash.f90 | 31 +- src/io/cli.f90 | 20 +- src/linalg.f90 | 175 ++++++++- src/new/envelopes.f90 | 46 +-- src/new/mod_inj_envelopes.f90 | 343 ++++++++++++------ src/optim.f90 | 4 +- src/phase_equilibria/mod_phase_equilibria.f90 | 75 ++-- .../mod_phase_equilibria.fypp | 15 +- tools/plot_pt.gnu | 14 +- tools/plot_px.gnu | 21 +- 12 files changed, 699 insertions(+), 265 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index 76f2dbc..c37337a 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -17,30 +17,40 @@ program main type(envelope) :: pt_bub, pt_dew, pt_hpl !! Shared 2ph-PT envelopes type(PTEnvel3), allocatable :: pt_bub_3(:), pt_dew_3(:) !! Shared 3ph-PT envelopes - type(injelope), allocatable :: px_bub(:), px_dew(:), px_hpl !! Shared 2ph-Px envelopes + type(injelope), allocatable :: px_bub(:), px_dew(:) !! Shared 2ph-Px envelopes + type(injelope) :: px_hpl(1) + + integer :: cli_error + logical :: run_px, run_3ph - ! real(pr) :: alpha=0.95 real(pr) :: alpha=0.0 real(pr) :: pt_bub_t0 = 180 real(pr) :: pt_dew_t0 = 180 ! Setup everything call setup - + + call cli%get(alpha, switch="--alpha0", error=cli_error) call get_z(alpha, z) + print *, "ALPHA: ", alpha + + call cli%get(run_3ph, switch="--three_phases", error=cli_error) + ! PT Envelopes call cpu_time(st) call pt_envelopes call cpu_time(et) print *, "PT: ", (et - st)*1000, "cpu ms" - - call exit + call cli%get(run_px, switch="--injection", error=cli_error) + + if (run_px) then ! PX Envelopes call cpu_time(st) call px_envelopes call cpu_time(et) print *, "PX: ", (et - st)*1000, "cpu ms" + end if contains subroutine setup !! Setup system @@ -83,7 +93,7 @@ subroutine pt_envelopes real(pr) :: t, p ! Temperature and pressure real(pr), allocatable :: k(:) ! K factors - integer :: n_points, icri(4), ncri, i + integer :: n_points, icri(4), ncri, i, i_max type(point), allocatable :: intersections(:), self_intersections(:) character(len=:), allocatable :: pt_case @@ -100,7 +110,8 @@ subroutine pt_envelopes ! Bubble envel ! ------------------------------------------------------------------------ ! call k_wilson_bubble(z, t_0=pt_bub_t0, p_end=0.5_pr, t=t, p=p, k=k) - bubble = bubble_temperature(z, p=1.0_pr, t0=300.0_pr) + print *, "Bubble PT" + bubble = bubble_temperature(z, p=1.0_pr, t0=pt_bub_t0) k = bubble%y/z t = bubble%t p = bubble%p @@ -115,6 +126,7 @@ subroutine pt_envelopes ! ======================================================================== ! Dew/AOP envelopes ! ------------------------------------------------------------------------ + print *, "Dew PT" t = pt_dew_t0 p = p_wilson(z, t) do while (p > 0.1) @@ -124,7 +136,7 @@ subroutine pt_envelopes k = 1/k_wilson(t, p) - dew = dew_temperature(z, p=0.1_pr, t0=300.0_pr) + dew = dew_temperature(z, p=0.00001_pr, t0=pt_dew_t0) k = dew%y/z t = dew%t p = dew%p @@ -139,18 +151,16 @@ subroutine pt_envelopes ! ======================================================================== ! HPLL Envelope ! ------------------------------------------------------------------------ - ! t = 700.0_pr - ! p = maxval([pt_bub%p, pt_dew%p])*1.5_pr - p = 300 - t = 550 ! pt_dew%t(minloc(pt_dew%p, dim=1))+50 + i_max = maxloc([pt_dew%p], dim=1) + + p = pt_dew%p(i_max) + 50 + t = pt_dew%t(i_max) + 50 - dew = hpl_temperature(z, p=p, t0=t) + dew = hpl_temperature(z, p=p, t0=t, y0=exp(pt_dew%logk(i_max, :)) * z) k = dew%y/z t = dew%t p = dew%p - print *, dew%iters, t, p - ! call find_hpl(t, p, k) call envelope2( & 3, nc, z, T, P, k, & @@ -169,6 +179,7 @@ subroutine pt_envelopes print *, style_bold // pt_case // style_reset ! ======================================================================== + if (run_3ph) then three_phase: block use envelopes, only: pt_three_phase_from_intersection @@ -183,13 +194,14 @@ subroutine pt_envelopes isolated: block use legacy_ar_models, only: nc, termo, z use envelopes, only: pt_envelope_three_phase + integer :: ncomp=1 real(pr) :: p, v, t real(pr) :: lnphix(nc), lnphiy(nc), & phase_x(nc), phase_y(nc), beta, x(2*nc+3), lnKx(nc), lnKy(nc) - beta = z(nc) + beta = z(ncomp) phase_y = 0 - phase_y(nc) = 1 + phase_y(ncomp) = 1 p = pt_bub%p(1) t = pt_bub%t(1) @@ -258,6 +270,7 @@ subroutine pt_envelopes ) end select end block three_phase + end if end subroutine subroutine px_envelopes @@ -271,9 +284,9 @@ subroutine px_envelopes px_hpl_line use envelopes, only: envelope, k_wilson, p_wilson use linalg, only: interpol, point, intersection - type(point), allocatable :: inter_dew_bub(:), self_inter_dew(:), self_inter_bub(:) - real(pr) :: t_tol = 2 + real(pr) :: t_tol = 2, p_hpl, a_hpl, z_hpl(nc), y_hpl(nc) + integer :: max_bub_p print *, style_underline // "----------" // style_reset print *, style_underline // "Px Regions" // style_reset @@ -283,85 +296,133 @@ subroutine px_envelopes ! Two phase envelopes ! ------------------------------------------------------------------------ print *, red // "Running Bubble" // style_reset - px_bub = px_two_phase_from_pt(t_inj, pt_bub, t_tol=5.0_pr) + px_bub = px_two_phase_from_pt(t_inj, pt_bub, alpha0=alpha, t_tol=5.0_pr) print *, blue // "Running Dew" // style_reset - px_dew = px_two_phase_from_pt(t_inj, pt_dew, t_tol=5.0_pr) + px_dew = px_two_phase_from_pt(t_inj, pt_dew, alpha0=alpha, t_tol=5.0_pr) print *, blue // "Running HPLL" // style_reset - px_hpl = px_hpl_line(0.99_pr, 1000.0_pr) + ! TODO: This is a dirty setup that barely works and should be fixed + + ! if (allocated(px_dew)) then + ! if (size(px_dew(1)%alpha) > 5) then + ! max_bub_p = maxloc(px_dew(1)%p, dim=1) + ! a_hpl = px_dew(1)%alpha(max_bub_p)*0.9 + ! p_hpl = px_dew(1)%p(max_bub_p) * 1.1 + ! call get_z(a_hpl, z_hpl) + ! y_hpl = exp(px_dew(1)%logk(max_bub_p+2, :)) * z + ! print *, a_hpl, p_hpl + ! px_hpl(1) = px_hpl_line(a_hpl, p_hpl, y0=y_hpl) + ! endif + ! if (size(px_bub(1)%alpha) > 5) then + ! max_bub_p = maxloc(px_bub(1)%p, dim=1) + ! a_hpl = px_bub(1)%alpha(max_bub_p)*1.1 + ! p_hpl = px_bub(1)%p(max_bub_p) + ! call get_z(a_hpl, z_hpl) + ! y_hpl = exp(px_bub(1)%logk(max_bub_p+2, :)) * z + ! px_hpl(1) = px_hpl_line(a_hpl, p_hpl, y0=y_hpl) + ! end if + ! end if + ! ======================================================================== ! ======================================================================== ! Three phase regions ! ------------------------------------------------------------------------ + if (run_3ph) then three_phase: block use dsp_lines, only: injelope, dsp_line_from_dsp_px - integer :: i, j, k type(injelope) :: px_bub_3, px_dew_3, px_branch_3(2) - type(injelope):: dsps(2) ! ===================================================================== ! Intersections between lines ! --------------------------------------------------------------------- - do i=1,size(px_dew) - do j=1,size(px_bub) - ! Go through each possible pair of envelopes to find DSPs - inter_dew_bub = intersection(& - px_dew(i)%alpha, px_dew(i)%p, & - px_bub(j)%alpha, px_bub(j)%p & - ) - do k=1,size(inter_dew_bub) - ! For each DSP found, calculate the two possible DSP lines - ! and the two three-phase branchs. - print *, "Intersection: ", inter_dew_bub(k) - dsps = dsp_line_from_dsp_px(& - inter_dew_bub(k), px_dew(i), px_bub(j) & - ) - px_branch_3 = px_three_phase_from_inter(& - inter_dew_bub(k), px_dew(i), px_bub(j) & - ) - end do - end do - end do - ! ===================================================================== + call calc_all_dsps(px_dew, px_bub, px_branch_3) + call calc_all_dsps(px_bub, px_hpl, px_branch_3) + call calc_all_dsps(px_dew, px_hpl, px_branch_3) + call calc_all_self_dsp(px_dew, px_branch_3) + call calc_all_self_dsp(px_bub, px_branch_3) ! ===================================================================== - ! Self intersections loops - ! --------------------------------------------------------------------- - do i=1,size(px_bub) - print *, "Checking Bub slef-intersections: ", i - self_inter_bub = intersection(px_bub(i)%alpha, px_bub(i)%p) - if (size(self_inter_bub) > 0) then - do j=1,size(self_inter_bub) - dsps = dsp_line_from_dsp_px(self_inter_bub(j), px_bub(i), px_bub(i)) - px_branch_3 = px_three_phase_from_inter(& - self_inter_bub(j), px_bub(i), px_bub(i) & - ) - end do - end if - end do - - do i=1,size(px_dew) - print *, "Checking Dew self-intersections: ", i - self_inter_dew = intersection(px_dew(i)%alpha, px_dew(i)%p) - if (size(self_inter_dew) > 0) then - do j=1,size(self_inter_dew) - dsps = dsp_line_from_dsp_px(self_inter_dew(j), px_dew(i), px_dew(i)) - px_branch_3 = px_three_phase_from_inter(& - self_inter_dew(j), px_dew(i), px_dew(i) & - ) - end do - end if - end do - ! ===================================================================== - + ! Isolated lines coming from PT lines print *, "Isolated Bub" - px_bub_3 = px_three_phase_from_pt(t_inj, pt_bub_3, t_tol) - print *, "Isolated Dew" - px_dew_3 = px_three_phase_from_pt(t_inj, pt_dew_3, t_tol) + if (allocated(pt_bub_3)) then + px_bub_3 = px_three_phase_from_pt(t_inj, pt_bub_3, 2*t_tol) + end if + + ! print *, "Isolated Dew" + if (allocated(pt_dew_3)) then + px_dew_3 = px_three_phase_from_pt(t_inj, pt_dew_3, t_tol) + end if end block three_phase + end if ! ======================================================================== end subroutine + + subroutine calc_all_dsps(px_1, px_2, px_out) + !! Calculate three-phase lines from DSPs between two lines. + !! + !! Find all the intersection points between two phase envelopes and + !! calculate all the pairs of three phase lines (\(P \alpha\) and DSP + !! lines) + use inj_envelopes, only: injelope, px_three_phase_from_inter + use dsp_lines, only: dsp_line_from_dsp_px + use linalg, only: point, intersection + type(injelope) :: px_1(:), px_2(:) + type(injelope) :: px_out(2) + type(injelope) :: dsps(2) + type(point), allocatable :: inter_1_2(:) + integer :: i, j, k + + do i=1,size(px_1) + do j=1,size(px_2) + ! Go through each possible pair of envelopes to find DSPs + inter_1_2 = intersection(& + px_1(i)%alpha, px_1(i)%p, & + px_2(j)%alpha, px_2(j)%p & + ) + do k=1,size(inter_1_2) + ! For each DSP found, calculate the two possible DSP lines + ! and the two three-phase branchs. + print *, "Intersection: ", inter_1_2(k) + dsps = dsp_line_from_dsp_px(& + inter_1_2(k), px_1(i), px_2(j) & + ) + px_out = px_three_phase_from_inter(& + inter_1_2(k), px_1(i), px_2(j) & + ) + end do + end do + end do + end subroutine + + subroutine calc_all_self_dsp(px, px_out) + !! Calculate three-phase lines from DSPs on a single line. + !! + !! From a single envelope find all self-crossing points an calculate + !! the corresponding three-phase lines for each one of them. Including + !! both Px as well as DSP lines. + use inj_envelopes, only: injelope, px_three_phase_from_inter + use dsp_lines, only: dsp_line_from_dsp_px + use linalg, only: point, intersection + type(injelope) :: px(:) + + type(injelope) :: px_out(2) + type(injelope) :: dsps(2) + type(point), allocatable :: inter(:) + integer :: i, j + + do i=1,size(px) + inter = intersection(px(i)%alpha, px(i)%p) + if (size(inter) > 0) then + do j=1,size(inter) + dsps = dsp_line_from_dsp_px(inter(j), px(i), px(i)) + px_out = px_three_phase_from_inter(& + inter(j), px(i), px(i) & + ) + end do + end if + end do + end subroutine end program main diff --git a/fpm.toml b/fpm.toml index 564609e..29a9abc 100644 --- a/fpm.toml +++ b/fpm.toml @@ -6,7 +6,7 @@ maintainer = "federico.benelli@mi.unc.edu.ar" copyright = "Copyright 2023, Federico E. Benelli" [build] -link = ["blas"] +link = ["lapack", "blas"] auto-executables = true auto-tests = true auto-examples = true @@ -30,4 +30,5 @@ implicit-typing = false # default: false [preprocess] [preprocess.cpp] -[preprocess.fypp] \ No newline at end of file +[preprocess.fypp] +suffixes = ["fypp"] diff --git a/src/Flash.f90 b/src/Flash.f90 index c72f566..ede8239 100644 --- a/src/Flash.f90 +++ b/src/Flash.f90 @@ -1,21 +1,36 @@ module phase_equilibria use constants, only: pr use legacy_ar_models, only: zTVTERMO, termo, n => nc, omg => w, tc, pc + use saturation_points, only: EquilibriaState implicit none contains - subroutine flash_pt(z, p, t, x, y, beta, its) + type(EquilibriaState) function flash_tp(z, t, p) real(pr), intent(in) :: z(:) !! Feed phase molar fractions - real(pr), intent(in) :: p !! Pressure [bar] real(pr), intent(in) :: t !! Temperature [K] - real(pr), intent(out) :: x(size(z)) !! Phase X molar fractions - real(pr), intent(out) :: y(size(z)) !! Phase Y molar fractions - real(pr), intent(out) :: beta !! Molar Y fraction - real(pr), optional, intent(out) :: its !! Number of iterations + real(pr), intent(in) :: p !! Pressure [bar] - real(pr) :: rho_x, rho_y - end subroutine + real(pr), :: x(size(z)) !! Phase X molar fractions + real(pr), :: y(size(z)) !! Phase Y molar fractions + real(pr), :: beta !! Molar Y fraction + real(pr), :: its !! Number of iterations + real(pr) :: rho_x, rho_y, v + + call flash("PT", .TRUE. z, t, p, v, x, y, rho_x, rho_y, beta, its) + + flash_tp = EquilibriaState(& + iters=its, y=y, x=x, beta=beta, rho_x=rho_x, rho_y=rho_y, t=t, p=p & + ) + if (rho_x < rho_y) then + flash_tp%x = y + flash_tp%y = x + + flash_tp%rho_x = rho_y + flash_tp%rho_y = rho_x + end if + + end function subroutine flash(spec, FIRST, z, t, p, v, x, y, rho_x, rho_y, beta, iter) ! Flash specification, eos id and number of compounds in the system diff --git a/src/io/cli.f90 b/src/io/cli.f90 index 667b847..531c708 100644 --- a/src/io/cli.f90 +++ b/src/io/cli.f90 @@ -22,8 +22,26 @@ subroutine setup_cli_envelopes(cli) switch_ab="-px", & help="Trace Px lines", & error=cli_error, & - required=.false.) + required=.false., & + def=".false.") + + call cli%add( & + switch="--three_phases", & + switch_ab="-3ph", & + help="Trace three phase lines", & + error=cli_error, & + required=.false., & + def=".false.") + + call cli%add( & + switch="--alpha0", & + switch_ab="-a", & + help="Initial alpha0, used for PT envelopes", & + error=cli_error, & + required=.false., & + def="0") call cli%parse(error=cli_error) + print *, cli_error if (cli_error /= 0) stop end subroutine diff --git a/src/linalg.f90 b/src/linalg.f90 index 951ff17..06e5a35 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -57,7 +57,9 @@ function intersect_two_lines(l1_x, l1_y, l2_x, l2_y) result(intersections) real(pr) :: s, t integer :: i, j - real(pr) :: x, y, xold = 9999, yold = 9999 + real(pr) :: x, y, xold, yold + xold = 9999 + yold = 9999 allocate (intersections(0)) @@ -115,6 +117,11 @@ function intersect_one_line(lx, ly) result(intersections) y3 => ly(j), y4 => ly(j - 1)) call intersects(x1, x2, x3, x4, y1, y2, y3, y4, s, t) + ! print *, "===", i, j,"===" + ! print *, x1, y1, x2, y2 + ! print *, x3, y3, x4, y4 + ! print *, s, t + ! print *, "---------------" if (0 <= s .and. s <= 1 .and. 0 <= t .and. t <= 1) then x = s*(x2 - x1) + x1 y = s*(y2 - y1) + y1 @@ -265,3 +272,169 @@ subroutine fcn(n, x, fvec, fjac, ldfjac, iflag) end subroutine end subroutine end module linalg + + +! module polyfit_mod +! +! implicit none +! interface +! subroutine sgelss (m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, info) +! integer :: m +! integer :: n +! integer :: nrhs +! real, dimension( lda, * ) :: a +! integer :: lda +! real, dimension( ldb, * ) :: b +! integer :: ldb +! real, dimension( * ) :: s +! real :: rcond +! integer :: rank +! real, dimension( * ) :: work +! integer :: lwork +! integer :: info +! end subroutine +! end interface +! +! contains +! +! function vander(x,N,increasing) result(v) +! real, intent(in) :: x(:) +! integer, intent(in) :: N +! logical, intent(in), optional :: increasing +! real :: v(size(x),N) +! +! integer :: i +! logical :: increasing_ +! +! increasing_ = .false. +! if (present(increasing)) increasing_ = increasing +! +! if (increasing_) then +! if (N > 0) v(:,1) = 1 +! if (N > 1) then +! v(:,2:) = spread(x,dim=2,ncopies=N-1) +! do i = 2, N-1 +! v(:,i+1) = v(:,i)*v(:,i+1) +! end do +! end if +! else +! if (N > 0) v(:,N) = 1 +! if (N > 1) then +! v(:,1:N-1) = spread(x,dim=2,ncopies=N-1) +! do i = N-1,2,-1 +! v(:,i-1) = v(:,i)*v(:,i-1) +! end do +! end if +! end if +! end function +! +! function polyfit(x,y,deg,rcond,rank,singular_values) result(p) +! real, intent(in) :: x(:), y(:) +! integer, intent(in) :: deg +! real, intent(in), optional :: rcond +! integer, intent(out), optional :: rank +! real, intent(out), allocatable, optional :: singular_values(:) +! real :: p(deg + 1) +! +! integer :: order +! real :: rcond_ +! real, allocatable :: lhs(:,:), rhs(:), scl(:),s(:), work(:) +! integer :: m, n, rank_, info, lwork +! +! order = deg + 1 +! +! ! check parameters +! if (deg < 0) error stop "expected deg >= 0" +! if (size(x) == 0) error stop "expected non-empty vector for x" +! if (size(x) /= size(y)) error stop "expected x and y to have same length" +! +! ! set rcond +! rcond_ = size(x)*epsilon(x) +! if (present(rcond)) rcond_ = rcond +! +! m = size(x) +! n = order +! allocate(lhs(m,n),scl(m)) +! +! ! set up least squares equation for powers of x +! lhs = vander(x,order) +! +! call print_mat(lhs) +! ! scale lhs to improve the condition number +! scl = sqrt(sum(lhs*lhs,dim=1)) +! print *, "scl = " +! call print_mat(spread(scl,dim=1,ncopies=m)) +! lhs = lhs/spread(scl,dim=1,ncopies=m) +! write(*,*) +! call print_mat(lhs) +! +! allocate(rhs,source=y) +! +! ! +! ! solve lstsq problem +! ! +! allocate(s(min(m,n))) +! +! ! perform workspace query +! lwork = -1 +! allocate(work(1)) +! call sgelss(m,n,1,lhs,m,rhs,m,s,rcond_,rank_,work,lwork,info) +! print *, "info = ", info +! +! ! resize workspace +! lwork = int(work(1)) +! print *, "lwork = ", lwork, work(1) +! deallocate(work) +! allocate(work(lwork)) +! +! print *, rhs +! ! compute minimum norm solutiong using SVD +! call sgelss(m,n,1,lhs,m,rhs,m,s,rcond_,rank_,work,lwork,info) +! if (info /= 0) then +! if (info < 0) write(*,*) 'sgelss: the ',abs(info),'-th argument had an illegal value' +! if (info > 0) write(*,*) 'sgelss: the algorithm for computing the SVD failed to converge' +! error stop +! end if +! call print_mat(lhs) +! print *, rhs +! +! if (rank_ /= order) then +! write(*,*) "polyfit may be poorly conditioned" +! end if +! +! ! broadcast scale coefficients +! p = rhs/scl +! +! if (present(rank)) rank = rank_ +! if (present(singular_values)) then +! if (allocated(singular_values)) deallocate(singular_values) +! allocate(singular_values,source=s) +! end if +! end function +! +! subroutine print_mat(v) +! real, intent(in) :: v(:,:) +! integer :: i, j +! do i = 1, size(v,1) +! print *, (v(i,j),j=1,size(v,2)) +! end do +! end subroutine +! +! end module +! +! program main +! +! use polyfit_mod +! implicit none +! +! real :: x(4), y(4) +! real :: p(3) +! +! x = [1.,2.,3.,4.] +! y = x**2 + 2*x - 4 +! +! p = polyfit(x,y,2) +! print *, p +! +! end program +! \ No newline at end of file diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index 8db7498..779ff99 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -102,7 +102,7 @@ subroutine update_specification(iter, passingcri, X, dF, ns, S, delS, dXdS) ! Selection of (the most changing) variable to be specified for the next point nsold = ns - ns = maxloc(abs(dXdS), dim=1) + ns = maxloc(abs(dXdS(:nc)), dim=1) if (maxval(abs(X(:nc))) < 0.2) then ns = maxloc(abs(dXdS(:nc)), dim=1) ! T and P not allowed to be chosen close to a critical point @@ -369,6 +369,8 @@ subroutine envelope2(ichoice, n, z, T, P, KFACT, & ! This will probably always e integer :: funit_output character(len=254) :: fname_env + real(pr) :: specifications(max_points, 3) + ! ======================================================================== ! OUTPUT file ! ------------------------------------------------------------------------ @@ -447,8 +449,12 @@ subroutine envelope2(ichoice, n, z, T, P, KFACT, & ! This will probably always e i=0 + open(69, file="specs" // str(env_number) // ".dat") do while (run) i = i + 1 + specifications(i, 1) = ns + specifications(i, 2) = S + specifications(i, 3) = delS if (i > max_points - 50) then exit end if @@ -556,42 +562,22 @@ subroutine envelope2(ichoice, n, z, T, P, KFACT, & ! This will probably always e ! to find a possible critical point (check if all lnK values ! change sign). extrapolation: block + ! use polyfit_mod, only: polyfit integer :: loc integer :: its real(pr) :: delta - real(pr) :: m(size(X)) + real(pr) :: m real(pr) :: max_lnK, max_lnK2, delta_lnK real(pr) :: delta_X(size(x)) real(pr) :: Xin(size(X)) its = 0 delta = delS - ! Variation of lnK based on deltaS - ! m = 37.0_pr * (X - Xold2)/(delta) - ! lnK_extrapolated = (delta) * m(:n) + X(:n) - lnK_extrapolated = X(:n) + 7 * delS * dXdS(:n) - ! print *, X(:n) - ! print *, lnK_extrapolated - - if (all((X(:n) * lnK_extrapolated < 0), dim=1)) then - print *, "Extrapol CP" - ! All lnK changed sign, so a CP is inminent - ! aproach it enough to force the jumping algorithm - do while( & - maxval(abs(X(:n))) >= 0.02 & - .and. all(X(:n)*lnK_extrapolated > 0, dim=1)& - ) - print *, its, "Getting to critical", & - exp(X(n+1)), exp(X(n+2)), maxval(abs(X(:n))), & - all(X(:n)*lnK_extrapolated > 0, dim=1) - its = its + 1 - delta_x = delta * m - X = delta_x + X - S = X(ns) - if (its > 10) exit - end do - passingcri = .true. + m = 2 + lnK_extrapolated = X(:n) + m * delS * dXdS(:n) + if (all(X(:n)*lnK_extrapolated < 0)) then + delS = 1.5 * m * delS end if end block extrapolation end if @@ -648,6 +634,10 @@ subroutine envelope2(ichoice, n, z, T, P, KFACT, & ! This will probably always e end if end do !----------------------------------------------------------- + write(69, *) specifications(:, 1) + write(69, *) specifications(:, 2) + write(69, *) specifications(:, 3) + close(69) n_points = i @@ -908,7 +898,7 @@ subroutine pt_envelope_three_phase(X0, spec_number, del_S0, envel) real(pr) :: Xnew(size(X0)) real(pr) :: dP, dT - del_S = sign(1.0_pr, del_S) * minval([ & + del_S = sign(1.5_pr, del_S) * minval([ & max(abs(sqrt(X(ns))/10), 0.1_pr), & abs(del_S)*3/iters & ] & diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index 41060b1..402906a 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -205,6 +205,9 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) end if XS(point, :) = X + call get_z(X(n+2), z) + write (funit_output, *) "SOL", iters, ns, X(n + 2), exp(X(n + 1)), & + X(:n), z call update_spec(X, ns, del_S, dF, dXdS) call fix_step_two_phases(X, ns, S, iters, del_S, dXdS) @@ -216,29 +219,32 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) integer :: max_changing fact = critical_fact ! 2.5 - Xnew = X + fact*dXdS*del_S - - K = X(:n) - Knew = Xnew(:n) + if (maxval(abs(X(:n))) < 1e-2) then + else + Xnew = X + fact*dXdS*del_S + K = X(:n) + Knew = Xnew(:n) - if (all(K*Knew < 0)) then - max_changing = maxloc(abs(K - Knew), dim=1) + if (all(K*Knew < 0)) then + print *, "CP", env_number + max_changing = maxloc(abs(K - Knew), dim=1) - dS_c = ( & - -K(max_changing)*(Xnew(max_changing) - X(max_changing)) & - /(Knew(max_changing) - K(max_changing)) & - ) + dS_c = ( & + -K(max_changing)*(Xnew(max_changing) - X(max_changing)) & + /(Knew(max_changing) - K(max_changing)) & + ) - Xnew = X + dXdS * dS_c - alpha_c = Xnew(n + 2) - pc = exp(Xnew(n + 1)) - cps = [cps, critical_point(t, pc, alpha_c)] - - del_S = dS_c + critical_multiplier * dS_c + Xnew = X + dXdS * dS_c + alpha_c = Xnew(n + 2) + pc = exp(Xnew(n + 1)) + cps = [cps, critical_point(t, pc, alpha_c)] + + del_S = dS_c + critical_multiplier * dS_c - write (funit_output, *) "" - write (funit_output, *) "" - end if + write (funit_output, *) "" + write (funit_output, *) "" + end if + endif end block detect_critical X = X + dXdS*del_S @@ -248,10 +254,6 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) print *, "Breaking: ", break_conditions(X, ns, S, del_S) exit enveloop end if - - call get_z(X(n+2), z) - write (funit_output, *) "SOL", iters, ns, X(n + 2), exp(X(n + 1)), & - X(:n), z end do enveloop write (funit_output, *) "" @@ -268,6 +270,8 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) close (funit_output) point = point - 1 + allocate(envels%t(point)) + envels%t = t envels%z = z_0 envels%z_inj = z_injection envels%logk = XS(:point, :n) @@ -362,8 +366,9 @@ function break_conditions(X, ns, S, del_S) p = exp(X(n + 1)) alpha = X(n + 2) - break_conditions = [ .false. & - ! p < 1e-15 .or. p > 5000, & + break_conditions = [ & + ! p > 50000, & + all(X(:n) < 1e-10) & ! abs(del_S) < 1e-3 & ] end function @@ -602,14 +607,14 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) exp(X(2*n + 1)), X(2*n + 3), X(:2*n), z XS(point, :) = X - call update_spec(X, ns, del_S, dF, dXdS) + call update_spec_three_phases(X, ns, del_S, dF, dXdS) fix_step: block real(pr) :: Xnew(size(X0)) real(pr) :: dP, dalpha del_S = sign(del_S_multiplier_three_phase, del_S)*minval([ & - max(abs(sqrt(X(ns))/10), 0.1_pr), & + max(abs(sqrt(X(ns))/5), 0.1_pr), & abs(del_S)*3/iters & ] & ) @@ -634,26 +639,30 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) Xnew(size(X0)), fact real(pr) :: pc, alpha_c, dS_c, dXdS_in(size(X0)) real(pr) :: bu_solve_tol - integer :: max_changing, i + integer :: max_changing, i_k logical :: found_critical fact = critical_fact found_critical = .false. - loop: do i = 0, 1 + Xnew = X + dXdS*del_S + do while (maxval(abs(Xnew(:2*n))) < 0.1) + print *, "increasing step" + del_S = 2*del_S + Xnew = X + dXdS*del_S + exit detect_critical + end do + + loop: do i_k = 0, 1 Xnew = X + fact*dXdS*del_S - K = X(i*n + 1:(i + 1)*n) - Knew = Xnew(i*n + 1:(i + 1)*n) + K = X(i_k*n + 1:(i_k + 1)*n) + Knew = Xnew(i_k*n + 1:(i_k + 1)*n) max_changing = maxloc(abs(Knew - K), dim=1) if (all(K*Knew < 0)) then - - print *, "" - print *, "CP ENV", env_number - print *, "" - + print *, "CRITICAL!" dS_c = ( & -k(max_changing) * (Xnew(ns) - X(ns)) & /(Knew(max_changing) - K(max_changing)) & @@ -664,7 +673,7 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) pc = exp(Xnew(2*n + 1)) cps = [cps, critical_point(t, pc, alpha_c)] - del_S = dS_c + 1.7_pr * dS_c + del_S = fact * dS_c write (funit_output, *) "" write (funit_output, *) "" @@ -724,10 +733,43 @@ function break_conditions_three_phases(X, ns, S, del_S) alpha = X(2*n + 2) break_conditions_three_phases = [ & - ! p < 0.001_pr .or. p > 5000, & - abs(del_S) < 1e-5 & + p > 50000 & + ! abs(del_S) < 1e-5 & ] end function + + subroutine update_spec_three_phases(X, ns, del_S, dF, dXdS) + real(pr), intent(in) :: X(:) + integer, intent(in out) :: ns + real(pr), intent(in out) :: del_S + real(pr), intent(in) :: dF(size(X), size(X)) + real(pr), intent(in out) :: dXdS(size(X)) + + real(pr) :: dFdS(size(X)) + integer :: ns_new + integer :: nc + + dFdS = 0 + dFdS(size(dFdS)) = 1 + dXdS = solve_system(dF, dFdS) + + nc = (size(X) - 3)/2 + + ns_new = maxloc(abs(dXdS(:2*nc)), dim=1) + ! if (any(abs(X(:2*nc)) < 0.1)) then + ! ns_new = maxloc(abs(dXdS(2:nc)), dim=1) + ! else + ! ns_new = maxloc(abs(dXdS), dim=1) + ! endif + + if (ns_new /= ns) then + ! translation of delS and dXdS to the new specification variable + del_S = dXdS(ns_new)*del_S + dXdS = dXdS/dXdS(ns_new) + ns = ns_new + end if + end subroutine + ! =========================================================================== subroutine get_z(alpha, z, dzda) @@ -783,7 +825,7 @@ function same_line(env1, env2) ! Initialization procedures ! --------------------------------------------------------------------------- function px_two_phase_from_pt(& - t_inj, pt_env_2, t_tol, del_S0) result(px_envels) + t_inj, pt_env_2, t_tol, alpha0, del_S0) result(px_envels) !! Calculate two phase Px envelopes at a given injection temperature. !! !! Given an injection temperature `t_inj` and a base PT envelope @@ -798,6 +840,7 @@ function px_two_phase_from_pt(& type(envelope), intent(in) :: pt_env_2 !! Base PT envelope real(pr), intent(in) :: t_tol !! Absolute temperature tolerance real(pr), optional, intent(in) :: del_S0 !! First point \(\Delta S\) + real(pr), optional, intent(in) :: alpha0 !! First point \(\alpha\) type(injelope), allocatable :: px_envels(:) !! Output Px envelope real(pr), allocatable :: ts_envel(:) !! Temperatures under tolerance @@ -813,6 +856,7 @@ function px_two_phase_from_pt(& type(injelope) :: px_envel allocate(px_envels(0)) + alpha = optval(alpha0, 0.0_pr) del_S = optval(del_S0, 0.1_pr) pold = 1e9 @@ -824,7 +868,7 @@ function px_two_phase_from_pt(& pt_env_2%p(idx), pt_env_2%p(idx + 1), & t_inj) - if (abs(p - pold) < 2) cycle + if (abs(p - pold) < 0.001) cycle pold = p k = exp(interpol( & @@ -832,7 +876,6 @@ function px_two_phase_from_pt(& pt_env_2%logk(idx, :), pt_env_2%logk(idx + 1, :), & t_inj)) - alpha = 0 X = [log(K), log(P), alpha] ns = size(X) @@ -872,9 +915,8 @@ function px_three_phase_from_pt(t_inj, pt_env_3, t_tol, del_S0) result(envel) del_S = optval(del_S0, 0.05_pr) pold = 0 - + do i_envel = 1, size(pt_env_3) - print *, i_envel associate(pt => pt_env_3(i_envel)) ts_envel = pack(pt%t, mask=abs(pt%t - t_inj) < t_tol) do i = 1, size(ts_envel) @@ -904,7 +946,7 @@ function px_three_phase_from_pt(t_inj, pt_env_3, t_tol, del_S0) result(envel) X = [log(Kx), log(Ky), log(P), alpha, beta] ns = size(X) - 1 - print *, "Running isolated PX" + print *, "Running isolated PX", alpha, P call injection_envelope_three_phase(X, ns, del_S, envel) end do end associate @@ -981,13 +1023,16 @@ function px_three_phase_from_inter(& call injection_envelope_three_phase(X, ns, del_S, envels(2)) end function - function px_hpl_line(alpha_0, p) - ! Find a HPLL PX line at a given pressure, starting from a given alpha + function px_hpl_line(alpha_0, p, y0) + !! Find a HPLL PX line at a given pressure, starting from a given alpha use legacy_ar_models, only: nc, termo use linalg, only: solve_system - real(pr), intent(in) :: alpha_0 !! Staring \(\alpha\) to search - real(pr), intent(in) :: p !! Pressure of HPLL + use saturation_points, only: EquilibriaState + real(pr), intent(in out) :: alpha_0 !! Staring \(\alpha\) to search + real(pr), intent(in out) :: p !! Pressure of HPLL + real(pr), optional :: y0(:) type(injelope) :: px_hpl_line !! Resulting HPLL line + type(injelope) :: px_hpl_line_1, px_hpl_line_2 !! intermediary HPLL lines real(pr) :: diff real(pr) :: lnfug_z(nc), lnfug_y(nc), & @@ -998,71 +1043,159 @@ function px_hpl_line(alpha_0, p) real(pr), allocatable :: x(:) real(pr) :: del_S0 - integer :: ns, ncomp - - find_hpl: do ncomp = nc, 1, -1 - alpha_in = alpha_0 - p_in = p - - y = 0 - y(ncomp) = 1.0 + integer :: ns, ncomp, npoints - diff = foo([alpha_in, p_in, y]) - return - do while (abs(diff) > 0.001 .and. p_in < 5000) - p_in = p_in + 10.0_pr - diff = foo([alpha_in, p_in, y]) - end do + integer :: i - if (p_in >= 5000) then - do while (abs(diff) > 0.001 .and. alpha_in > 0) - alpha_in = alpha_in - 0.01_pr - diff = foo([alpha_in, p_in, y]) - end do - end if + type(EquilibriaState) :: hpl_state - if (.true.) then!(alpha_in > 0 .or. p_in < 5000) then - call optim - print *, alpha_in, p_in - k = 1/exp(lnfug_y - lnfug_z) + do i=0,int(p),50 + hpl_state = hpl_alpha_pressure(z, t, p-i, alpha_0, y0) + call get_z(alpha_0, z) + ! print *, sum(z*hpl_state%y/hpl_state%x)- 1, alpha_0, hpl_state%p + print *, hpl_state%x(:5) + print *, hpl_state%y(:5) + if (maxval(abs(hpl_state%x - hpl_state%y)) > 0.1) exit + end do + + x = hpl_state%x + y = hpl_state%y + print *, hpl_state%iters + print *, x(:5) + print *, y(:5) + + X = [log(x/y), log(hpl_state%p), alpha_0] + del_S0 = 0.5_pr + ns = size(X) - 1 + call injection_envelope(X, ns, del_S0, px_hpl_line_1) + + ! if (size(px_hpl_line_1%alpha) > 1) then + ! ! If the first line give results, calculate the other segment + ! ! going to the other direction + ! X = [px_hpl_line_1%logk(1, :), log(px_hpl_line_1%p(1)), px_hpl_line_1%alpha(1)] + ! del_S0 = -0.5_pr + ! call injection_envelope(X, ns, del_S0, px_hpl_line_2) + + ! px_hpl_line%alpha = [px_hpl_line_2%alpha, px_hpl_line_1%alpha] + ! px_hpl_line%t = [px_hpl_line_2%t, px_hpl_line_1%t] + ! px_hpl_line%p = [px_hpl_line_2%p, px_hpl_line_1%p] + + ! npoints = size(px_hpl_line_1%alpha) + size(px_hpl_line_2%alpha) + ! allocate(px_hpl_line%logk(npoints, nc)) + ! px_hpl_line%logk(:size(px_hpl_line_2%alpha), :) = px_hpl_line_2%logk + ! px_hpl_line%logk(size(px_hpl_line_2%alpha)+1:, :) = px_hpl_line_1%logk + ! else + ! px_hpl_line = px_hpl_line_1 + ! end if + end function - X = [log(K), log(P_in), alpha_in] + type(EquilibriaState) function hpl_alpha_pressure(n, t, p0, a0, y0, max_inner_its) + !! HPLL \(\alpha\) and P calculation on a specified T. + use saturation_points, only: EquilibriaState + use envelopes, only: k_wilson + use stdlib_optval, only: optval + use legacy_ar_models, only: nc, termo + real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] + real(pr), intent(in) :: t !! Temperature [K] + real(pr), intent(in) :: p0 !! Pressure [bar] + real(pr), intent(in out) :: a0 !! \(\alpha_0\) + real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations + + real(pr) :: vy, vz + + real(pr) :: k(size(n)), y(size(n)), z(size(n)), dzda(size(n)), lnk(size(n)) + real(pr) :: lnfug_y(size(n)), dlnphi_dz(size(n), size(n)) + real(pr) :: lnfug_z(size(n)), dlnphi_dy(size(n), size(n)) + real(pr) :: dkda(size(n)), dydz(size(n)) + + real(pr) :: x(2), f, df(2), step + + integer :: i, its, inner_its + integer :: max_iterations = 100 + real(pr) :: step_tol=0.1_pr, tol=1.e-5_pr + real(pr) :: p + ! ======================================================================= + ! Handle arguments + ! ----------------------------------------------------------------------- + if (size(n) /= nc) call exit(1) + z = n/sum(n) + inner_its = optval(inner_its, 50) + + ! Initiliaze with K-wilson factors + if (present(y0)) then + y = y0 + k = y/z + else + k = k_wilson(t, p0) + y = k * z + end if + ! ======================================================================= + + ! ======================================================================= + ! Solve point + ! ----------------------------------------------------------------------- + x = [a0, p0] + block + use optimization, only: nm_opt + integer :: stat + real(pr) :: step(2) + x = [a0, p] + step = [0.01_pr, 1._pr] + call nm_opt(foo, x, stat, step) + end block + + ! do its=1,max_iterations + ! f = foo(x) + ! df = dfoo(x) + ! x = x - 5*df + ! if (maxval(abs(df)) < tol) exit + ! end do + + a0 = x(1) + p = x(2) + + call get_z(a0, z, dzda) + call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, fugn=dlnphi_dz) + call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, fugn=dlnphi_dy) + print *, vy, vz + hpl_alpha_pressure = EquilibriaState(its, y, z, t, p) + ! ======================================================================= + contains + function foo(x) + real(pr) :: x(:) + real(pr) :: alpha + real(pr) :: pressure + real(pr) :: foo + + alpha = x(1) + pressure = x(2) + call get_z(alpha, z, dzda) + call termo(nc, 1, 4, t, pressure, y, vy, philog=lnfug_y, fugn=dlnphi_dz) + call termo(nc, 1, 4, t, pressure, z, vz, philog=lnfug_z, fugn=dlnphi_dy) + + lnk = lnfug_z - lnfug_y + k = exp(lnk) + foo = (sum(z*k) - 1)**2 + end function - del_S0 = 0.1_pr - ns = size(X) - 1 - call injection_envelope(X, ns, del_S0, px_hpl_line) - exit find_hpl - end if - end do find_hpl + function dfoo(x) + real(pr) :: x(2) + real(pr) :: dfoo(2) - contains - function foo(x) result(f) - real(pr) :: x(:) + real(pr) :: dx(2) + real(pr) :: fdx real(pr) :: f - if (x(1) > 1) x(1) = 0.97_pr - ! alpha_in = x(1) - p_in = x(2) - y = abs(x(3:)) + dx = 0 + dx(1) = 0.001_pr + + dfoo(1) = (foo(x+dx) - foo(x))/dx(1) - call get_z(alpha_in, z) - call termo(nc, 4, 1, t, p_in, z, v, philog=lnfug_z, dlphip=dlnphi_dp_z, fugn=dlnphi_dz) - call termo(nc, 4, 1, t, p_in, y, v, philog=lnfug_y, dlphip=dlnphi_dp_y, fugn=dlnphi_dy) - f = abs((log(z(ncomp)) + lnfug_z(ncomp)) - (log(y(ncomp)) + lnfug_y(ncomp))) + dx = 0 + dx(2) = 0.01_pr + dfoo(2) = (foo(x+dx) - foo(x))/dx(2) end function - - subroutine optim - use optimization, only: nm_opt - real(pr) :: x0(size(z) + 2) - integer :: stat - real(pr) :: step(size(x0)) - - x0 = [alpha_in, p_in, y] - step = 0.1 - step(1) = 0.1 - step(2) = -100 - call nm_opt(foo, x0, stat, step) - end subroutine end function ! =========================================================================== -end module \ No newline at end of file +end module diff --git a/src/optim.f90 b/src/optim.f90 index 3c3861e..83599b7 100644 --- a/src/optim.f90 +++ b/src/optim.f90 @@ -108,9 +108,9 @@ subroutine nelmin(fn, n, start, xmin, ynewlo, reqmin, step, konvge, kcount, & real(kind=rk), parameter :: ecoeff = 2.0D+00 real(kind=rk), parameter :: eps = 0.001D+00 interface - function fn(x) + function fn(xx) import rk - real(rk) :: x(:) + real(rk) :: xx(:) real(rk) :: fn end function end interface diff --git a/src/phase_equilibria/mod_phase_equilibria.f90 b/src/phase_equilibria/mod_phase_equilibria.f90 index a0453cf..7fbc40c 100644 --- a/src/phase_equilibria/mod_phase_equilibria.f90 +++ b/src/phase_equilibria/mod_phase_equilibria.f90 @@ -16,14 +16,17 @@ module saturation_points integer :: max_iterations = 100 real(pr) :: step_tol = 0.1_pr contains - type(EquilibriaState) function hpl_pressure(n, t, p0, y0) + + type(EquilibriaState) function hpl_pressure(n, t, p0, y0, max_inner_its) !! Hpl pressure calculation function. !! !! Calculates the hpl pressure of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: t !! Temperature [K] real(pr), optional, intent(in) :: p0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: p, vy, vz @@ -33,7 +36,7 @@ type(EquilibriaState) function hpl_pressure(n, t, p0, y0) real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -41,14 +44,15 @@ type(EquilibriaState) function hpl_pressure(n, t, p0, y0) z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (p0)) p = p0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 k = y/z else k = k_wilson(t, p) - y = k * z + y = k * z end if ! ======================================================================= @@ -58,8 +62,10 @@ type(EquilibriaState) function hpl_pressure(n, t, p0, y0) do its=1,max_iterations call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphip=dlnphi_dp_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphip=dlnphi_dp_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 p = p/2.0_pr call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphip=dlnphi_dp_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphip=dlnphi_dp_z) @@ -73,19 +79,22 @@ type(EquilibriaState) function hpl_pressure(n, t, p0, y0) end do p = p - step y = z*k - if (abs(step) < tol) exit + if (abs(step) < tol .and. its > 3) exit end do hpl_pressure = EquilibriaState(its, y, z, t, p) ! ======================================================================= end function - type(EquilibriaState) function hpl_temperature(n, p, t0, y0) + + type(EquilibriaState) function hpl_temperature(n, p, t0, y0, max_inner_its) !! Hpl temperature calculation function. !! !! Calculates the hpl temperature of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: p !! Temperature [K] real(pr), optional, intent(in) :: t0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: t, vy, vz @@ -95,7 +104,7 @@ type(EquilibriaState) function hpl_temperature(n, p, t0, y0) real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -103,7 +112,8 @@ type(EquilibriaState) function hpl_temperature(n, p, t0, y0) z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (t0)) t = t0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 @@ -120,8 +130,10 @@ type(EquilibriaState) function hpl_temperature(n, p, t0, y0) do its=1,max_iterations call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphit=dlnphi_dt_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphit=dlnphi_dt_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 t = t - 5.0_pr call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphit=dlnphi_dt_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphit=dlnphi_dt_z) @@ -141,14 +153,16 @@ type(EquilibriaState) function hpl_temperature(n, p, t0, y0) ! ======================================================================= end function - type(EquilibriaState) function bubble_pressure(n, t, p0, y0) + type(EquilibriaState) function bubble_pressure(n, t, p0, y0, max_inner_its) !! Bubble pressure calculation function. !! !! Calculates the bubble pressure of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: t !! Temperature [K] real(pr), optional, intent(in) :: p0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: p, vy, vz @@ -158,7 +172,7 @@ type(EquilibriaState) function bubble_pressure(n, t, p0, y0) real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -166,7 +180,8 @@ type(EquilibriaState) function bubble_pressure(n, t, p0, y0) z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (p0)) p = p0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 @@ -183,8 +198,10 @@ type(EquilibriaState) function bubble_pressure(n, t, p0, y0) do its=1,max_iterations call termo(nc, -1, 4, t, p, y, vy, philog=lnfug_y, dlphip=dlnphi_dp_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphip=dlnphi_dp_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 p = p/2.0_pr call termo(nc, -1, 4, t, p, y, vy, philog=lnfug_y, dlphip=dlnphi_dp_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphip=dlnphi_dp_z) @@ -203,14 +220,17 @@ type(EquilibriaState) function bubble_pressure(n, t, p0, y0) bubble_pressure = EquilibriaState(its, y, z, t, p) ! ======================================================================= end function - type(EquilibriaState) function bubble_temperature(n, p, t0, y0) + + type(EquilibriaState) function bubble_temperature(n, p, t0, y0, max_inner_its) !! Bubble temperature calculation function. !! !! Calculates the bubble temperature of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: p !! Temperature [K] real(pr), optional, intent(in) :: t0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: t, vy, vz @@ -220,7 +240,7 @@ type(EquilibriaState) function bubble_temperature(n, p, t0, y0) real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -228,7 +248,8 @@ type(EquilibriaState) function bubble_temperature(n, p, t0, y0) z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (t0)) t = t0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 @@ -245,8 +266,10 @@ type(EquilibriaState) function bubble_temperature(n, p, t0, y0) do its=1,max_iterations call termo(nc, -1, 4, t, p, y, vy, philog=lnfug_y, dlphit=dlnphi_dt_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphit=dlnphi_dt_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 t = t - 5.0_pr call termo(nc, -1, 4, t, p, y, vy, philog=lnfug_y, dlphit=dlnphi_dt_y) call termo(nc, 1, 4, t, p, z, vz, philog=lnfug_z, dlphit=dlnphi_dt_z) @@ -266,14 +289,16 @@ type(EquilibriaState) function bubble_temperature(n, p, t0, y0) ! ======================================================================= end function - type(EquilibriaState) function dew_pressure(n, t, p0, y0) + type(EquilibriaState) function dew_pressure(n, t, p0, y0, max_inner_its) !! Dew pressure calculation function. !! !! Calculates the dew pressure of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: t !! Temperature [K] real(pr), optional, intent(in) :: p0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: p, vy, vz @@ -283,7 +308,7 @@ type(EquilibriaState) function dew_pressure(n, t, p0, y0) real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -291,7 +316,8 @@ type(EquilibriaState) function dew_pressure(n, t, p0, y0) z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (p0)) p = p0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 @@ -308,8 +334,10 @@ type(EquilibriaState) function dew_pressure(n, t, p0, y0) do its=1,max_iterations call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphip=dlnphi_dp_y) call termo(nc, -1, 4, t, p, z, vz, philog=lnfug_z, dlphip=dlnphi_dp_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 p = p/2.0_pr call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphip=dlnphi_dp_y) call termo(nc, -1, 4, t, p, z, vz, philog=lnfug_z, dlphip=dlnphi_dp_z) @@ -328,14 +356,17 @@ type(EquilibriaState) function dew_pressure(n, t, p0, y0) dew_pressure = EquilibriaState(its, y, z, t, p) ! ======================================================================= end function - type(EquilibriaState) function dew_temperature(n, p, t0, y0) + + type(EquilibriaState) function dew_temperature(n, p, t0, y0, max_inner_its) !! Dew temperature calculation function. !! !! Calculates the dew temperature of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: p !! Temperature [K] real(pr), optional, intent(in) :: t0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: t, vy, vz @@ -345,7 +376,7 @@ type(EquilibriaState) function dew_temperature(n, p, t0, y0) real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -353,7 +384,8 @@ type(EquilibriaState) function dew_temperature(n, p, t0, y0) z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (t0)) t = t0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 @@ -370,8 +402,10 @@ type(EquilibriaState) function dew_temperature(n, p, t0, y0) do its=1,max_iterations call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphit=dlnphi_dt_y) call termo(nc, -1, 4, t, p, z, vz, philog=lnfug_z, dlphit=dlnphi_dt_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 t = t - 5.0_pr call termo(nc, 1, 4, t, p, y, vy, philog=lnfug_y, dlphit=dlnphi_dt_y) call termo(nc, -1, 4, t, p, z, vz, philog=lnfug_z, dlphit=dlnphi_dt_z) @@ -390,5 +424,4 @@ type(EquilibriaState) function dew_temperature(n, p, t0, y0) dew_temperature = EquilibriaState(its, y, z, t, p) ! ======================================================================= end function - end module diff --git a/src/phase_equilibria/mod_phase_equilibria.fypp b/src/phase_equilibria/mod_phase_equilibria.fypp index a19e18f..adf976d 100644 --- a/src/phase_equilibria/mod_phase_equilibria.fypp +++ b/src/phase_equilibria/mod_phase_equilibria.fypp @@ -23,7 +23,7 @@ contains #: if val=="pressure" #: set var="p" #: set spec="t" - #: else + #: elif val=="temperature" #: set var="t" #: set spec="p" #: endif @@ -37,14 +37,16 @@ contains #: set ROOTy=1 #: set ROOTz=1 #: endif - type(EquilibriaState) function ${inci}$_${val}$(n, ${spec}$, ${var}$0, y0) + type(EquilibriaState) function ${inci}$_${val}$(n, ${spec}$, ${var}$0, y0, max_inner_its) !! ${inci.title()}$ ${val}$ calculation function. !! !! Calculates the ${inci}$ ${val}$ of a multicomponent mixture. + use stdlib_optval, only: optval real(pr), intent(in) :: n(:) !! Composition vector [moles / molar fraction] real(pr), intent(in) :: ${spec}$ !! Temperature [K] real(pr), optional, intent(in) :: ${var}$0 !! Initial pressure [bar] real(pr), optional, intent(in) :: y0(:) !! Initial composition + integer, optional, intent(in) :: max_inner_its(:) !! Inner iterations real(pr) :: ${var}$, vy, vz @@ -54,7 +56,7 @@ contains real(pr) :: f, step - integer :: its + integer :: its, inner_its ! ======================================================================= ! Handle arguments @@ -62,7 +64,8 @@ contains z = n/sum(n) if (size(n) /= nc) call exit(1) if (present (${var}$0)) ${var}$ = ${var}$0 - + inner_its = optval(inner_its, 50) + ! Initiliaze with K-wilson factors if (present(y0)) then y = y0 @@ -83,11 +86,13 @@ contains do its=1,max_iterations call termo(nc, ${ROOTy}$, 4, t, p, y, vy, philog=lnfug_y, dlphi${var}$=dlnphi_d${var}$_y) call termo(nc, ${ROOTz}$, 4, t, p, z, vz, philog=lnfug_z, dlphi${var}$=dlnphi_d${var}$_z) + inner_its = 0 do while (any(isnan(lnfug_y)) .and. t > 0) + inner_its = inner_its + 1 #:if var == "p" ${var}$ = ${var}$/2.0_pr - #:else + #:elif var == "t" ${var}$ = ${var}$ - 5.0_pr #:endif call termo(nc, ${ROOTy}$, 4, t, p, y, vy, philog=lnfug_y, dlphi${var}$=dlnphi_d${var}$_y) diff --git a/tools/plot_pt.gnu b/tools/plot_pt.gnu index fedc70d..27d4fb2 100644 --- a/tools/plot_pt.gnu +++ b/tools/plot_pt.gnu @@ -10,13 +10,13 @@ set mytics set xlabel "Temperature [K]" set ylabel "Pressure [bar]" -plot "fenvelopes_output/env-2ph-PT_1.dat" u 4:5 ,\ - "fenvelopes_output/env-2ph-PT_2.dat" u 4:5 ,\ - "fenvelopes_output/env-2ph-PT_3.dat" u 4:5 ,\ - "fenvelopes_output/env-3ph-PT_4.dat" u 4:5 ,\ - "fenvelopes_output/env-3ph-PT_5.dat" u 4:5 ,\ - "fenvelopes_output/env-3ph-PT_6.dat" u 4:5 ,\ - "fenvelopes_output/env-3ph-PT_7.dat" u 4:5 ,\ +plot "fenvelopes_output/env-2ph-PT_1.dat" u 4:5 w lp,\ + "fenvelopes_output/env-2ph-PT_2.dat" u 4:5 w lp,\ + "fenvelopes_output/env-2ph-PT_3.dat" u 4:5 w lp,\ + "fenvelopes_output/env-3ph-PT_4.dat" u 4:5 w lp,\ + "fenvelopes_output/env-3ph-PT_5.dat" u 4:5 w lp,\ + "fenvelopes_output/env-3ph-PT_6.dat" u 4:5 w lp,\ + "fenvelopes_output/env-3ph-PT_7.dat" u 4:5 w lp,\ #plot "fenvelopes_output/env-2ph-PT_1.dat" u 4:5 w l lc "black" ,\ # "fenvelopes_output/env-2ph-PT_2.dat" u 4:5 w l lc "blue" ,\ diff --git a/tools/plot_px.gnu b/tools/plot_px.gnu index d73011b..ab0f90c 100644 --- a/tools/plot_px.gnu +++ b/tools/plot_px.gnu @@ -33,13 +33,18 @@ set style line liq_3ph lc rgb "red" dt 7 lw 1.5 pxs_2 = system("ls -d ./fenvelopes_output/env-2ph-PX* | xargs") pxs_3 = system("ls -d ./fenvelopes_output/env-2ph-PX* | xargs") -plot "fenvelopes_output/env-2ph-PX_1.dat" u 4:5 w lp ls bub_2ph t "2ph-bub",\ - "fenvelopes_output/env-2ph-PX_2.dat" u 4:5 w lp ls dew_2ph t "2ph-dew", \ - "fenvelopes_output/env-2ph-PX_3.dat" u 4:5 w lp ls dew_2ph t "2ph-dew", \ - "fenvelopes_output/env-2ph-PX_4.dat" u 4:5 w lp ls dew_2ph t "2ph-dew", \ - "fenvelopes_output/env-3ph-PX_3.dat" u 4:5 w lp ls bub_3ph t "3ph-bub", \ - "fenvelopes_output/env-3ph-PX_4.dat" u 4:5 w lp ls dew_3ph t "3ph-dew", \ - "fenvelopes_output/env-3ph-PX_5.dat" u 4:5 w lp ls bub_3ph t "3ph-bub", \ - "fenvelopes_output/env-3ph-PX_6.dat" u 4:5 w lp ls dew_3ph t "3ph-dew", \ +plot "fenvelopes_output/env-2ph-PX_1.dat" u 4:5 w l ls bub_2ph t "2ph-bub",\ + "fenvelopes_output/env-2ph-PX_1.dat" index "critical" u 1:2 w p ls bub_2ph t "",\ + "fenvelopes_output/env-2ph-PX_2.dat" u 4:5 w l ls dew_2ph t "2ph-dew", \ + "fenvelopes_output/env-2ph-PX_2.dat" index "critical" u 1:2 w p ls bub_2ph t "",\ + "fenvelopes_output/env-2ph-PX_3.dat" u 4:5 w l ls dew_2ph t "2ph-dew", \ + "fenvelopes_output/env-2ph-PX_4.dat" u 4:5 w l ls dew_2ph t "2ph-dew", \ + "fenvelopes_output/env-3ph-PX_3.dat" u 4:5 w l ls bub_3ph t "3ph-bub", \ + "fenvelopes_output/env-3ph-PX_4.dat" u 4:5 w l ls dew_3ph t "3ph-dew", \ + "fenvelopes_output/env-3ph-PX_5.dat" u 4:5 w l ls bub_3ph t "3ph-bub", \ + "fenvelopes_output/env-3ph-PX_6.dat" u 4:5 w l ls dew_3ph t "3ph-dew", \ + "fenvelopes_output/env-3ph-PX_7.dat" u 4:5 w l ls bub_3ph t "3ph-bub", \ + "fenvelopes_output/env-3ph-PX_8.dat" u 4:5 w l ls dew_3ph t "3ph-dew" + pause mouse close