From 150a755575ffcba2db6b031b4929c6e57dda4371 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:19:36 -0300 Subject: [PATCH 01/29] feat(DSP): Calculation of DSP lines --- src/new/mod_dsp_lines.f90 | 515 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 515 insertions(+) create mode 100644 src/new/mod_dsp_lines.f90 diff --git a/src/new/mod_dsp_lines.f90 b/src/new/mod_dsp_lines.f90 new file mode 100644 index 0000000..42697a8 --- /dev/null +++ b/src/new/mod_dsp_lines.f90 @@ -0,0 +1,515 @@ +module dsp_lines + !! Module to calculate DSP lines + use constants, only: pr, R + use dtypes, only: envelope, critical_point + use inj_envelopes, only: update_spec, injection_case, z_0, z_injection, & + get_z, injelope + use linalg, only: solve_system, interpol, full_newton + use progress_bar_module, only: progress_bar + + implicit none + ! =========================================================================== + ! Parameters + ! --------------------------------------------------------------------------- + integer :: env_number = 0 !! Number of calculated envelope + integer :: max_iters = 50 !! Maximum number of iterations for a newton step + integer, parameter :: max_points = 1000 !! Maximum number of points for each envelope + character(len=255) :: FE_LOG + + ! Two-phase settings + real(pr) :: del_S_multiplier = 2.0_pr + real(pr) :: max_dp=50.0_pr + real(pr) :: max_dalpha=0.01_pr + real(pr) :: critical_multiplier = 2.0_pr + + ! Three phase parameters + real(pr) :: del_S_multiplier_three_phase = 1.7_pr + real(pr) :: critical_fact = 3.0_pr + ! =========================================================================== +contains + subroutine from_nml(filepath) + use legacy_ar_models, only: nc + character(len=*), intent(in) :: filepath + integer :: funit + + namelist /nml_px/ T, z_0, z_injection, injection_case + + allocate (z_0(nc), z_injection(nc)) + + open (newunit=funit, file=filepath) + read (funit, nml=nml_px) + close (funit) + + z_injection = z_injection/sum(z_injection) + end subroutine + + ! =========================================================================== + ! Three-phases + ! --------------------------------------------------------------------------- + subroutine F_injection_three_phases(Xvars, ns, S, F, dF) + !! Function to solve at each point of a three phase envelope. + !! + !! The vector of variables X corresponds to: + !! \( X = [lnKx_i, lnKy_i lnP, \alpha, T] \) + !! + !! While the equations are: + !! + !! \( F = [ + !! lnKx_i - ln \phi_i(x, P, T) + ln \phi_i(w, P, T), + !! lnKy_i - ln \phi_i(y, P, T) + ln \phi_i(w, P, T), + !! \sum_{i=1}^N (w_i) - 1, + !! \sum_{i=1}^N (x_i - y_i), + !! X_{ns} - S + !! ] \) + use legacy_ar_models, only: TERMO + use iso_fortran_env, only: error_unit + real(pr), intent(in) :: Xvars(:) !! Vector of variables + integer, intent(in) :: ns !! Number of specification + real(pr), intent(in) :: S !! Specification value + real(pr), intent(out) :: F(size(Xvars)) !! Vector of functions valuated + real(pr), intent(out) :: df(size(Xvars), size(Xvars)) !! Jacobian matrix + + ! Xvars variables + real(pr) :: z((Size(Xvars)-3)/2) + real(pr) :: Kx((Size(Xvars)-3)/2) + real(pr) :: Ky((Size(Xvars)-3)/2) + real(pr) :: P + real(pr) :: alpha + real(pr) :: T + + real(pr) :: beta=1 + + ! Main phase 1 variables + real(pr) :: Vx + real(pr), dimension((Size(Xvars)-3)/2) :: x, lnfug_x, dlnphi_dt_x, dlnphi_dp_x + real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_x + + ! Main phase 2 variables + real(pr) :: Vy + real(pr), dimension((Size(Xvars)-3)/2) :: y, lnfug_y, dlnphi_dt_y, dlnphi_dp_y + real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_y + + ! Incipient phase variables + real(pr) :: Vw + real(pr), dimension((Size(Xvars)-3)/2) :: w, lnfug_w, dlnphi_dt_w, dlnphi_dp_w + real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_w + + ! Derivative of z wrt alpha + real(pr) :: dzda((Size(Xvars)-3)/2), dwda((Size(Xvars)-3)/2) + + ! Derivative of w wrt beta + real(pr) :: dwdb((Size(Xvars)-3)/2) + + real(pr) :: dwdKx((Size(Xvars)-3)/2), dxdKx((Size(Xvars)-3)/2), dydKx((Size(Xvars)-3)/2) + real(pr) :: dwdKy((Size(Xvars)-3)/2), dxdKy((Size(Xvars)-3)/2), dydKy((Size(Xvars)-3)/2) + + integer :: i, j, n + + n = (Size(Xvars)-3)/2 + + ! Setting variables + Kx = exp(Xvars(1:n)) + Ky = exp(Xvars(n + 1:2*n)) + P = exp(Xvars(2*n + 1)) + alpha = Xvars(2*n + 2) + T = Xvars(2*n + 3) + + call get_z(alpha, z, dzda) + + w = z/(beta*Ky + (1 - beta)*Kx) + x = w*Kx + y = w*Ky + + call TERMO( & + n, 0, 4, T, P, x, Vx, lnfug_x, dlnphi_dp_x, dlnphi_dt_x, dlnphi_dn_x & + ) + call TERMO( & + n, 0, 4, T, P, y, Vy, lnfug_y, dlnphi_dp_y, dlnphi_dt_y, dlnphi_dn_y & + ) + call TERMO( & + n, 0, 4, T, P, w, Vw, lnfug_w, dlnphi_dp_w, dlnphi_dt_w, dlnphi_dn_w & + ) + + F(1:n) = Xvars(1:n) + lnfug_x - lnfug_w + F(n + 1:2*n) = Xvars(n + 1:2*n) + lnfug_y - lnfug_w + + F(2*n + 1) = sum(w) - 1 + F(2*n + 2) = sum(x - y) + F(2*n + 3) = Xvars(ns) - S + + df = 0 + dwda = 1.0_pr/(beta*Ky + (1 - beta)*Kx)*dzda + dwdb = z*(Kx - Ky)/((1 - beta)*Kx + beta*Ky)**2 + + dwdKx = -z*(1 - beta)/(Ky*beta + (1 - beta)*Kx)**2 + dxdKx = Kx*dwdKx + w + dydKx = Ky*dwdKx + + dwdKy = -z*(beta)/(Ky*beta + (1 - beta)*Kx)**2 + dxdKy = Kx*dwdKy + dydKy = Ky*dwdKy + w + + do i = 1, n + do j = 1, n + df(i, j) = Kx(j)*(dlnphi_dn_x(i, j)*dxdKx(j) & + - dlnphi_dn_w(i, j)*dwdKx(j)) + df(i + n, j) = Kx(j)*(dlnphi_dn_y(i, j)*dydKx(j) & + - dlnphi_dn_w(i, j)*dwdKx(j)) + + df(i, j + n) = Ky(j)*(dlnphi_dn_x(i, j)*dxdKy(j) & + - dlnphi_dn_w(i, j)*dwdKy(j)) + df(i + n, j + n) = Ky(j)*(dlnphi_dn_y(i, j)*dydKy(j) & + - dlnphi_dn_w(i, j)*dwdKy(j)) + end do + + df(i, i) = df(i, i) + 1 + df(i + n, i + n) = df(i + n, i + n) + 1 + df(i, 2*n + 2) = sum( & + Kx*dlnphi_dn_x(i, :)*dwda - dlnphi_dn_w(i, :)*dwda & + ) + df(i + n, 2*n + 2) = sum(Ky*dlnphi_dn_y(i, :)*dwda & + - dlnphi_dn_w(i, :)*dwda) + df(i, 2*n + 3) = sum(Kx*dlnphi_dn_x(i, :)*dwdb & + - dlnphi_dn_w(i, :)*dwdb) + df(i + n, 2*n + 3) = sum(Ky*dlnphi_dn_y(i, :)*dwdb & + - dlnphi_dn_w(i, :)*dwdb) + df(2*n + 1, i) = Kx(i)*dwdKx(i) + df(2*n + 1, i + n) = Ky(i)*dwdKy(i) + + df(2*n + 2, i) = Kx(i)*dxdKx(i) - Kx(i)*dydKx(i) + df(2*n + 2, i + n) = Ky(i)*dxdKy(i) - Ky(i)*dydKy(i) + end do + + ! Derivatives wrt P + df(:n, 2*n + 1) = P*(dlnphi_dp_x - dlnphi_dp_w) + df(n + 1:2*n, 2*n + 1) = P*(dlnphi_dp_y - dlnphi_dp_w) + + ! Derivatives wrt alpha + df(2*n + 1, 2*n + 2) = sum(dwda) + df(2*n + 2, 2*n + 2) = sum(Kx*dwda - Ky*dwda) + + ! Derivatives wrt T + df(:n, 2*n + 3) = T*(dlnphi_dp_x - dlnphi_dp_w) + df(n + 1:2*n, 2*n + 1) = T*(dlnphi_dp_y - dlnphi_dp_w) + + ! Derivatives wrt Xs + df(2*n + 3, :) = 0 + df(2*n + 3, ns) = 1 + end subroutine + + subroutine dsp_line(X0, spec_number, del_S0, envels) + use constants, only: ouput_path + use io, only: str + !! Subroutine to calculate Px phase envelopes via continuation method. + !! Three phases version. + real(pr), intent(in) :: X0(:) !! Vector of variables + integer, intent(in) :: spec_number !! Number of specification + real(pr), intent(in) :: del_S0 !! \(\Delta S_0\) + type(injelope), intent(out) :: envels !! Calculated envelopes + + type(critical_point), allocatable :: cps(:) + + real(pr) :: X(size(X0)) + integer :: ns + real(pr) :: S + real(pr) :: XS(max_points, size(X0)) + real(pr) :: del_S + + real(pr) :: F(size(X0)), dF(size(X0), size(X0)), dXdS(size(X0)) + + integer :: point, iters, n + integer :: i + integer :: funit_output + character(len=254) :: fname_env + + allocate (cps(0)) + X = X0 + n = (size(X0) - 3)/2 + ns = spec_number + S = X(ns) + del_S = del_S0 + + ! ====================================================================== + ! Output file + ! ---------------------------------------------------------------------- + env_number = env_number + 1 + + write (fname_env, *) env_number + fname_env = "env-3ph-DSP"//"_"//trim(adjustl(fname_env)) + fname_env = trim(adjustl(ouput_path))//trim(fname_env)//".dat" + + open (newunit=funit_output, file=fname_env) + write (funit_output, *) "#", T + write (funit_output, *) "STAT", " iters", " ns", " alpha", " P", & + " T", (" lnKx"//str(i), i=1,n), (" lnKy"//str(i), i=1,n) + write (funit_output, *) "X0", iters, ns, X(2*n + 2), exp(X(2*n + 1)), & + X(2*n + 3), X(:2*n) + ! ====================================================================== + + enveloop: do point = 1, max_points + call progress_bar(point, max_points, advance=.false.) + call full_newton(F_injection_three_phases, iters, X, ns, S, max_iters, F, dF) + if (iters >= max_iters) then + call progress_bar(point, max_points, advance=.true.) + print *, "Breaking: Above max iterations" + exit enveloop + end if + + write (funit_output, *) "SOL", iters, ns, X(2*n + 2), & + exp(X(2*n + 1)), X(2*n + 3), X(:2*n) + XS(point, :) = X + + call update_spec(X, ns, del_S, dF, dXdS) + + fix_step: block + end block fix_step + + detect_critical: block + real(pr) :: K((size(X0) - 3)/2), Knew((size(X0) - 3)/2), & + Xnew(size(X0)), fact + real(pr) :: pc, alpha_c, dS_c, dXdS_in(size(X0)) + integer :: max_changing, i + fact = critical_fact + + loop: do i = 0, 1 + Xnew = X + fact*dXdS*del_S + + K = X(i*n + 1:(i + 1)*n) + Knew = Xnew(i*n + 1:(i + 1)*n) + + max_changing = maxloc(abs(Knew - K), dim=1) + + if (all(K*Knew < 0)) then + dS_c = ( & + -k(max_changing)*(Xnew(ns) - X(ns)) & + /(Knew(max_changing) - K(max_changing)) & + ) + + Xnew = X + dXdS*dS_c + alpha_c = Xnew(2*n + 2) + pc = exp(Xnew(2*n + 1)) + cps = [cps, critical_point(t, pc, alpha_c)] + + ! del_S = dS_c + del_S ! * fact + ! del_S = del_S * fact + del_S = dS_c - sign(1.7_pr, dS_c)*dS_c + + write (funit_output, *) "" + write (funit_output, *) "" + exit loop + end if + end do loop + end block detect_critical + + if (x(2*n + 3) > 1 .or. (x(2*n+3) < 0)) then + call progress_bar(point, max_points, .true.) + print *, "Breaking: positive 𝜷" + exit enveloop + end if + + X = X + dXdS*del_S + S = X(ns) + if (any(break_conditions_three_phases(X, ns, S, del_S)) .and. point > 10) then + call progress_bar(point, max_points, .true.) + print *, "Breaking: ", break_conditions_three_phases(X, ns, S, del_S) + exit enveloop + end if + end do enveloop + + ! point = point - 1 + + write (funit_output, *) "" + write (funit_output, *) "" + write (funit_output, *) "#critical" + if (size(cps) > 0) then + do i = 1, size(cps) + write (funit_output, *) cps(i)%alpha, cps(i)%p + end do + else + write (funit_output, *) "NaN NaN" + end if + + close (funit_output) + envels%z = z_0 + envels%z_inj = z_injection + envels%logk = XS(:point, :n) + envels%alpha = XS(:point, n + 2) + envels%p = exp(XS(:point, n + 1)) + envels%critical_points = cps + end subroutine + + subroutine update_spec(X, ns, S, del_S, dXdS) + real(pr), intent(in) :: X(:) + integer, intent(in) :: ns + real(pr), intent(in) :: S + real(pr), intent(in) :: del_S + real(pr), intent(in out) :: dXdS(size(X)) + real(pr) :: Xnew(size(X)) + real(pr) :: dP, dT + + del_S = sign(del_S_multiplier_three_phase, del_S)*minval([ & + max(abs(X(ns)/10), 0.1_pr), & + abs(del_S)*3/iters & + ] & + ) + + if (injection_case == "dilution") del_S = 50*del_S + + Xnew = X + dXdS*del_S + dP = exp(Xnew(2*n + 1)) - exp(X(2*n + 1)) + dalpha = exp(Xnew(2*n + 3)) - exp(X(2*n + 3)) + + do while (abs(dP) > max_dp .or. abs(dalpha) > max_dalpha) + dXdS = dXdS/2.0_pr + + Xnew = X + dXdS*del_S + dP = exp(Xnew(2*n + 1)) - exp(X(2*n + 1)) + dalpha = Xnew(2*n + 2) - X(2*n + 2) + end do + end subroutine + + function break_conditions_dsp_line(X, ns, S, del_S) + !! Set of conditions to break the tracing. + real(pr) :: X(:) !! Variables vector + integer :: ns !! Number of specification + real(pr) :: S !! Value of specification + real(pr) :: del_S + + integer :: n + real(pr) :: p, alpha + logical, allocatable :: break_conditions_three_phases(:) + + n = (size(X) - 3)/2 + p = exp(X(2*n + 1)) + alpha = X(2*n + 2) + + break_conditions_three_phases = [ & + ! p < 0.001_pr .or. p > 5000, & + abs(del_S) < 1e-5 & + ] + end function + ! =========================================================================== + + subroutine get_z(alpha, z, dzda) + !! Calculate the fluid composition based on an amount of addition + !! of second fluid. + !! + !! The injection can be considered as two kinds of injection: + !! - Displacement: \( z = \alpha z_i + (1-\alpha) z_0 \) + !! - Addition: \( z = \frac{\alpha z_i + (1-\alpha) z_0}{\sum_{i=1}^N \alpha z_i + (1-\alpha) z_0} \) + real(pr), intent(in) :: alpha !! Addition percentaje \( \alpha \) + real(pr), intent(out) :: z(size(z_0)) !! New composition + real(pr), optional, intent(out) :: dzda(size(z_0)) !! Derivative wrt \(\alpha\) + + select case (injection_case) + case ("displace") + z = (z_injection*alpha + (1.0_pr - alpha)*z_0) + if (present(dzda)) dzda = z_injection - z_0 + case ("dilute") + z = (z_injection*alpha + z_0)/sum(z_injection*alpha + z_0) + if (present(dzda)) dzda = -(alpha*z_injection + z_0) & + *sum(z_injection)/sum(alpha*z_injection + z_0)**2 & + + z_injection/sum(alpha*z_injection + z_0) + case default + z = (z_injection*alpha + (1.0_pr - alpha)*z_0) + if (present(dzda)) dzda = z_injection - z_0 + end select + end subroutine + + function get_case(dew_envel, bub_envel) result(n_case) + type(injelope), intent(in) :: dew_envel + type(injelope), intent(in) :: bub_envel + integer :: n_case + end function + + function remove_duplicates(envels) result(clean_envels) + !! From a set of envelopes check if they are the same line + class(injelope) :: envels(:) + type(injelope), allocatable :: clean_envels(:) + + if (size(envels) /= 1) then + else + clean_envels = envels + end if + end function + + function same_line(env1, env2) + !! + class(injelope), intent(in) :: env1, env2 + logical :: same_line + end function + + ! =========================================================================== + ! Initialization procedures + ! --------------------------------------------------------------------------- + function px_three_phase_from_inter(& + inter, px_1, px_2, del_S0, beta0 & + ) result(envels) + use legacy_ar_models, only: nc + use stdlib_optval, only: optval + use linalg, only: point, interpol + type(point), intent(in) :: inter + type(injelope), intent(in) :: px_1, px_2 + type(injelope) :: envels(2) + real(pr), optional :: del_S0 + real(pr), optional :: beta0 + + integer :: i, j + + real(pr) :: lnKx(nc), lnKy(nc), alpha, p, beta, X(2*nc+3) + real(pr) :: phase_y(nc), phase_x(nc) + real(pr) :: del_S + real(pr) :: z(nc), dzda(nc) + integer :: ns + + del_S = optval(del_S0, -0.05_pr) + beta = optval(beta0, 0.99_pr) + + i = inter%i + j = inter%j + + alpha = inter%x + p = inter%y + + lnKx = interpol( & + px_1%alpha(i), px_1%alpha(i + 1), & + px_1%logk(i, :), px_1%logk(i + 1, :), & + alpha & + ) + + lnKy = interpol( & + px_2%alpha(j), px_2%alpha(j + 1), & + px_2%logk(j, :), px_2%logk(j + 1, :), & + alpha & + ) + + call get_z(alpha, z, dzda) + + ! Bubble line composition + phase_y = exp(lnKy)*z + ! Dew line composition + phase_x = exp(lnKx)*z + + ns = 2*nc + 3 + + ! ================================================================== + ! Line with incipient phase gas + ! ------------------------------------------------------------------ + print *, "Three Phase: Gas" + lnKx = log(phase_x/phase_y) + lnKy = log(z/phase_y) + X = [lnKx, lnKy, log(p), alpha, beta] + call injection_envelope_three_phase(X, ns, del_S, envels(1)) + ! ================================================================== + + ! ================================================================== + ! Line with incipient phase liquid + ! ------------------------------------------------------------------ + print *, "Three Phase: Liquid" + lnKx = log(phase_y/phase_x) + lnKy = log(z/phase_x) + X = [lnKx, lnKy, log(p), alpha, beta] + call injection_envelope_three_phase(X, ns, del_S, envels(2)) + end function + ! =========================================================================== +end module From 8f60a0b6f7d7196e7817fe3ec3d50a872f880756 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:19:59 -0300 Subject: [PATCH 02/29] cleaning of DSP mod --- src/new/mod_dsp_lines.f90 | 358 ++++++++++---------------------------- 1 file changed, 96 insertions(+), 262 deletions(-) diff --git a/src/new/mod_dsp_lines.f90 b/src/new/mod_dsp_lines.f90 index 42697a8..798d83f 100644 --- a/src/new/mod_dsp_lines.f90 +++ b/src/new/mod_dsp_lines.f90 @@ -12,41 +12,20 @@ module dsp_lines ! Parameters ! --------------------------------------------------------------------------- integer :: env_number = 0 !! Number of calculated envelope - integer :: max_iters = 50 !! Maximum number of iterations for a newton step + integer :: max_iters = 1000 !! Maximum number of iterations for a newton step integer, parameter :: max_points = 1000 !! Maximum number of points for each envelope character(len=255) :: FE_LOG - ! Two-phase settings - real(pr) :: del_S_multiplier = 2.0_pr - real(pr) :: max_dp=50.0_pr - real(pr) :: max_dalpha=0.01_pr - real(pr) :: critical_multiplier = 2.0_pr - - ! Three phase parameters - real(pr) :: del_S_multiplier_three_phase = 1.7_pr - real(pr) :: critical_fact = 3.0_pr + real(pr) :: del_S_multiplier = 1.0_pr + real(pr) :: max_da=0.1_pr + real(pr) :: max_dp=20.0_pr + real(pr) :: max_dT=10.0_pr ! =========================================================================== contains - subroutine from_nml(filepath) - use legacy_ar_models, only: nc - character(len=*), intent(in) :: filepath - integer :: funit - - namelist /nml_px/ T, z_0, z_injection, injection_case - - allocate (z_0(nc), z_injection(nc)) - - open (newunit=funit, file=filepath) - read (funit, nml=nml_px) - close (funit) - - z_injection = z_injection/sum(z_injection) - end subroutine - ! =========================================================================== ! Three-phases ! --------------------------------------------------------------------------- - subroutine F_injection_three_phases(Xvars, ns, S, F, dF) + subroutine dsp_line_F(Xvars, ns, S, F, dF) !! Function to solve at each point of a three phase envelope. !! !! The vector of variables X corresponds to: @@ -70,38 +49,32 @@ subroutine F_injection_three_phases(Xvars, ns, S, F, dF) real(pr), intent(out) :: df(size(Xvars), size(Xvars)) !! Jacobian matrix ! Xvars variables - real(pr) :: z((Size(Xvars)-3)/2) real(pr) :: Kx((Size(Xvars)-3)/2) real(pr) :: Ky((Size(Xvars)-3)/2) real(pr) :: P real(pr) :: alpha real(pr) :: T - real(pr) :: beta=1 - - ! Main phase 1 variables + ! Incipient phase 1 variables real(pr) :: Vx real(pr), dimension((Size(Xvars)-3)/2) :: x, lnfug_x, dlnphi_dt_x, dlnphi_dp_x real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_x - ! Main phase 2 variables + ! Incipient phase 2 variables real(pr) :: Vy real(pr), dimension((Size(Xvars)-3)/2) :: y, lnfug_y, dlnphi_dt_y, dlnphi_dp_y real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_y - ! Incipient phase variables - real(pr) :: Vw - real(pr), dimension((Size(Xvars)-3)/2) :: w, lnfug_w, dlnphi_dt_w, dlnphi_dp_w - real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_w + ! Main phase variables + real(pr) :: Vz + real(pr), dimension((Size(Xvars)-3)/2) :: z, lnfug_z, dlnphi_dt_z, dlnphi_dp_z + real(pr), dimension((Size(Xvars)-3)/2, (Size(Xvars)-3)/2) :: dlnphi_dn_z ! Derivative of z wrt alpha - real(pr) :: dzda((Size(Xvars)-3)/2), dwda((Size(Xvars)-3)/2) + real(pr) :: dzda((Size(Xvars)-3)/2) - ! Derivative of w wrt beta - real(pr) :: dwdb((Size(Xvars)-3)/2) - - real(pr) :: dwdKx((Size(Xvars)-3)/2), dxdKx((Size(Xvars)-3)/2), dydKx((Size(Xvars)-3)/2) - real(pr) :: dwdKy((Size(Xvars)-3)/2), dxdKy((Size(Xvars)-3)/2), dydKy((Size(Xvars)-3)/2) + real(pr) :: dzdKx((Size(Xvars)-3)/2), dxdKx((Size(Xvars)-3)/2), dydKx((Size(Xvars)-3)/2) + real(pr) :: dzdKy((Size(Xvars)-3)/2), dxdKy((Size(Xvars)-3)/2), dydKy((Size(Xvars)-3)/2) integer :: i, j, n @@ -110,91 +83,61 @@ subroutine F_injection_three_phases(Xvars, ns, S, F, dF) ! Setting variables Kx = exp(Xvars(1:n)) Ky = exp(Xvars(n + 1:2*n)) + P = exp(Xvars(2*n + 1)) alpha = Xvars(2*n + 2) - T = Xvars(2*n + 3) + T = exp(Xvars(2*n + 3)) call get_z(alpha, z, dzda) + ! if (any(z < 0)) return - w = z/(beta*Ky + (1 - beta)*Kx) - x = w*Kx - y = w*Ky + x = Kx * z + y = Ky * z call TERMO( & n, 0, 4, T, P, x, Vx, lnfug_x, dlnphi_dp_x, dlnphi_dt_x, dlnphi_dn_x & - ) + ) call TERMO( & n, 0, 4, T, P, y, Vy, lnfug_y, dlnphi_dp_y, dlnphi_dt_y, dlnphi_dn_y & - ) + ) call TERMO( & - n, 0, 4, T, P, w, Vw, lnfug_w, dlnphi_dp_w, dlnphi_dt_w, dlnphi_dn_w & - ) + n, 0, 4, T, P, z, Vz, lnfug_z, dlnphi_dp_z, dlnphi_dt_z, dlnphi_dn_z & + ) - F(1:n) = Xvars(1:n) + lnfug_x - lnfug_w - F(n + 1:2*n) = Xvars(n + 1:2*n) + lnfug_y - lnfug_w + F(1:n) = log(Kx) + lnfug_x - lnfug_z + F(n + 1:2*n) = log(Ky) + lnfug_y - lnfug_z - F(2*n + 1) = sum(w) - 1 - F(2*n + 2) = sum(x - y) + F(2*n + 1) = sum(x - z) + F(2*n + 2) = sum(y - z) F(2*n + 3) = Xvars(ns) - S - df = 0 - dwda = 1.0_pr/(beta*Ky + (1 - beta)*Kx)*dzda - dwdb = z*(Kx - Ky)/((1 - beta)*Kx + beta*Ky)**2 - - dwdKx = -z*(1 - beta)/(Ky*beta + (1 - beta)*Kx)**2 - dxdKx = Kx*dwdKx + w - dydKx = Ky*dwdKx - - dwdKy = -z*(beta)/(Ky*beta + (1 - beta)*Kx)**2 - dxdKy = Kx*dwdKy - dydKy = Ky*dwdKy + w - - do i = 1, n - do j = 1, n - df(i, j) = Kx(j)*(dlnphi_dn_x(i, j)*dxdKx(j) & - - dlnphi_dn_w(i, j)*dwdKx(j)) - df(i + n, j) = Kx(j)*(dlnphi_dn_y(i, j)*dydKx(j) & - - dlnphi_dn_w(i, j)*dwdKx(j)) - - df(i, j + n) = Ky(j)*(dlnphi_dn_x(i, j)*dxdKy(j) & - - dlnphi_dn_w(i, j)*dwdKy(j)) - df(i + n, j + n) = Ky(j)*(dlnphi_dn_y(i, j)*dydKy(j) & - - dlnphi_dn_w(i, j)*dwdKy(j)) - end do + do i=1,n + do j=1,n + df(i, j) = x(j) * dlnphi_dn_x(i, j) + df(i+n, j+n) = y(j) * dlnphi_dn_y(i, j) + end do df(i, i) = df(i, i) + 1 - df(i + n, i + n) = df(i + n, i + n) + 1 - df(i, 2*n + 2) = sum( & - Kx*dlnphi_dn_x(i, :)*dwda - dlnphi_dn_w(i, :)*dwda & - ) - df(i + n, 2*n + 2) = sum(Ky*dlnphi_dn_y(i, :)*dwda & - - dlnphi_dn_w(i, :)*dwda) - df(i, 2*n + 3) = sum(Kx*dlnphi_dn_x(i, :)*dwdb & - - dlnphi_dn_w(i, :)*dwdb) - df(i + n, 2*n + 3) = sum(Ky*dlnphi_dn_y(i, :)*dwdb & - - dlnphi_dn_w(i, :)*dwdb) - df(2*n + 1, i) = Kx(i)*dwdKx(i) - df(2*n + 1, i + n) = Ky(i)*dwdKy(i) - - df(2*n + 2, i) = Kx(i)*dxdKx(i) - Kx(i)*dydKx(i) - df(2*n + 2, i + n) = Ky(i)*dxdKy(i) - Ky(i)*dydKy(i) + df(i+n,i+n) = df(i+n, i+n) + 1 + + ! Derivatives wrt alpha + df(i, 2*n+2) = sum(dlnphi_dn_x(i, :) * Kx * dzda - dlnphi_dn_z(i, :) * dzda) + df(i+n, 2*n+2) = sum(dlnphi_dn_y(i, :) * Ky * dzda - dlnphi_dn_z(i, :) * dzda) end do - ! Derivatives wrt P - df(:n, 2*n + 1) = P*(dlnphi_dp_x - dlnphi_dp_w) - df(n + 1:2*n, 2*n + 1) = P*(dlnphi_dp_y - dlnphi_dp_w) + df(:n, 2*n+1) = P * (dlnphi_dp_x - dlnphi_dp_z) + df(n+1:2*n, 2*n+1) = P * (dlnphi_dp_y - dlnphi_dp_z) + + df(:n, 2*n+3) = T * (dlnphi_dt_x - dlnphi_dt_z) + df(n+1:2*n, 2*n+3) = T * (dlnphi_dt_y - dlnphi_dt_z) - ! Derivatives wrt alpha - df(2*n + 1, 2*n + 2) = sum(dwda) - df(2*n + 2, 2*n + 2) = sum(Kx*dwda - Ky*dwda) + df(2*n+1, :n) = x + df(2*n+2, n+1:2*n) = y - ! Derivatives wrt T - df(:n, 2*n + 3) = T*(dlnphi_dp_x - dlnphi_dp_w) - df(n + 1:2*n, 2*n + 1) = T*(dlnphi_dp_y - dlnphi_dp_w) + df(2*n+1, 2*n+2) = sum(Kx * dzda - dzda) + df(2*n+2, 2*n+2) = sum(Ky * dzda - dzda) - ! Derivatives wrt Xs - df(2*n + 3, :) = 0 - df(2*n + 3, ns) = 1 + df(2*n+3, ns) = 1 end subroutine subroutine dsp_line(X0, spec_number, del_S0, envels) @@ -239,85 +182,38 @@ subroutine dsp_line(X0, spec_number, del_S0, envels) fname_env = trim(adjustl(ouput_path))//trim(fname_env)//".dat" open (newunit=funit_output, file=fname_env) - write (funit_output, *) "#", T + write (funit_output, *) "#", env_number write (funit_output, *) "STAT", " iters", " ns", " alpha", " P", & " T", (" lnKx"//str(i), i=1,n), (" lnKy"//str(i), i=1,n) - write (funit_output, *) "X0", iters, ns, X(2*n + 2), exp(X(2*n + 1)), & - X(2*n + 3), X(:2*n) + write (funit_output, *) "X0", iters, ns, & + X(2*n + 2), exp(X(2*n + 1)), exp(X(2*n + 3)), X(:2*n) ! ====================================================================== enveloop: do point = 1, max_points call progress_bar(point, max_points, advance=.false.) - call full_newton(F_injection_three_phases, iters, X, ns, S, max_iters, F, dF) + call full_newton(dsp_line_F, iters, X, ns, S, max_iters, F, dF) if (iters >= max_iters) then call progress_bar(point, max_points, advance=.true.) print *, "Breaking: Above max iterations" exit enveloop end if - write (funit_output, *) "SOL", iters, ns, X(2*n + 2), & - exp(X(2*n + 1)), X(2*n + 3), X(:2*n) + write (funit_output, *) "SOL", iters, ns, & + X(2*n + 2), exp(X(2*n + 1)), exp(X(2*n + 3)), X(:2*n) XS(point, :) = X call update_spec(X, ns, del_S, dF, dXdS) - - fix_step: block - end block fix_step - - detect_critical: block - real(pr) :: K((size(X0) - 3)/2), Knew((size(X0) - 3)/2), & - Xnew(size(X0)), fact - real(pr) :: pc, alpha_c, dS_c, dXdS_in(size(X0)) - integer :: max_changing, i - fact = critical_fact - - loop: do i = 0, 1 - Xnew = X + fact*dXdS*del_S - - K = X(i*n + 1:(i + 1)*n) - Knew = Xnew(i*n + 1:(i + 1)*n) - - max_changing = maxloc(abs(Knew - K), dim=1) - - if (all(K*Knew < 0)) then - dS_c = ( & - -k(max_changing)*(Xnew(ns) - X(ns)) & - /(Knew(max_changing) - K(max_changing)) & - ) - - Xnew = X + dXdS*dS_c - alpha_c = Xnew(2*n + 2) - pc = exp(Xnew(2*n + 1)) - cps = [cps, critical_point(t, pc, alpha_c)] - - ! del_S = dS_c + del_S ! * fact - ! del_S = del_S * fact - del_S = dS_c - sign(1.7_pr, dS_c)*dS_c - - write (funit_output, *) "" - write (funit_output, *) "" - exit loop - end if - end do loop - end block detect_critical - - if (x(2*n + 3) > 1 .or. (x(2*n+3) < 0)) then - call progress_bar(point, max_points, .true.) - print *, "Breaking: positive 𝜷" - exit enveloop - end if + call fix_step(iters, X, ns, S, del_S, dXdS) X = X + dXdS*del_S S = X(ns) - if (any(break_conditions_three_phases(X, ns, S, del_S)) .and. point > 10) then + if (any(break_conditions_dsp_line(X, ns, S, del_S)) .and. point > 10) then call progress_bar(point, max_points, .true.) - print *, "Breaking: ", break_conditions_three_phases(X, ns, S, del_S) + print *, "Breaking: ", break_conditions_dsp_line(X, ns, S, del_S) exit enveloop end if end do enveloop - ! point = point - 1 - write (funit_output, *) "" write (funit_output, *) "" write (funit_output, *) "#critical" @@ -338,16 +234,20 @@ subroutine dsp_line(X0, spec_number, del_S0, envels) envels%critical_points = cps end subroutine - subroutine update_spec(X, ns, S, del_S, dXdS) + subroutine fix_step(iters, X, ns, S, del_S, dXdS) + integer, intent(in) :: iters real(pr), intent(in) :: X(:) integer, intent(in) :: ns real(pr), intent(in) :: S - real(pr), intent(in) :: del_S + real(pr), intent(in out) :: del_S real(pr), intent(in out) :: dXdS(size(X)) real(pr) :: Xnew(size(X)) - real(pr) :: dP, dT + real(pr) :: da, dP, dT + integer :: n + + n = (size(X) - 3)/2 - del_S = sign(del_S_multiplier_three_phase, del_S)*minval([ & + del_S = sign(del_S_multiplier, del_S)*minval([ & max(abs(X(ns)/10), 0.1_pr), & abs(del_S)*3/iters & ] & @@ -356,15 +256,16 @@ subroutine update_spec(X, ns, S, del_S, dXdS) if (injection_case == "dilution") del_S = 50*del_S Xnew = X + dXdS*del_S + da = (Xnew(2*n + 2)) - (X(2*n + 2)) dP = exp(Xnew(2*n + 1)) - exp(X(2*n + 1)) - dalpha = exp(Xnew(2*n + 3)) - exp(X(2*n + 3)) + dT = exp(Xnew(2*n + 3)) - exp(X(2*n + 3)) - do while (abs(dP) > max_dp .or. abs(dalpha) > max_dalpha) + do while (abs(da) > max_da .or. abs(dP) > max_dP .or. abs(dT) > max_dT) dXdS = dXdS/2.0_pr - Xnew = X + dXdS*del_S + da = (Xnew(2*n + 2)) - (X(2*n + 2)) dP = exp(Xnew(2*n + 1)) - exp(X(2*n + 1)) - dalpha = Xnew(2*n + 2) - X(2*n + 2) + dT = exp(Xnew(2*n + 3)) - exp(X(2*n + 3)) end do end subroutine @@ -377,139 +278,72 @@ function break_conditions_dsp_line(X, ns, S, del_S) integer :: n real(pr) :: p, alpha - logical, allocatable :: break_conditions_three_phases(:) + logical, allocatable :: break_conditions_dsp_line(:) n = (size(X) - 3)/2 p = exp(X(2*n + 1)) alpha = X(2*n + 2) - break_conditions_three_phases = [ & + break_conditions_dsp_line = [ .false.& ! p < 0.001_pr .or. p > 5000, & - abs(del_S) < 1e-5 & + ! abs(del_S) < 1e-5 & ] end function ! =========================================================================== - subroutine get_z(alpha, z, dzda) - !! Calculate the fluid composition based on an amount of addition - !! of second fluid. - !! - !! The injection can be considered as two kinds of injection: - !! - Displacement: \( z = \alpha z_i + (1-\alpha) z_0 \) - !! - Addition: \( z = \frac{\alpha z_i + (1-\alpha) z_0}{\sum_{i=1}^N \alpha z_i + (1-\alpha) z_0} \) - real(pr), intent(in) :: alpha !! Addition percentaje \( \alpha \) - real(pr), intent(out) :: z(size(z_0)) !! New composition - real(pr), optional, intent(out) :: dzda(size(z_0)) !! Derivative wrt \(\alpha\) - - select case (injection_case) - case ("displace") - z = (z_injection*alpha + (1.0_pr - alpha)*z_0) - if (present(dzda)) dzda = z_injection - z_0 - case ("dilute") - z = (z_injection*alpha + z_0)/sum(z_injection*alpha + z_0) - if (present(dzda)) dzda = -(alpha*z_injection + z_0) & - *sum(z_injection)/sum(alpha*z_injection + z_0)**2 & - + z_injection/sum(alpha*z_injection + z_0) - case default - z = (z_injection*alpha + (1.0_pr - alpha)*z_0) - if (present(dzda)) dzda = z_injection - z_0 - end select - end subroutine - - function get_case(dew_envel, bub_envel) result(n_case) - type(injelope), intent(in) :: dew_envel - type(injelope), intent(in) :: bub_envel - integer :: n_case - end function - - function remove_duplicates(envels) result(clean_envels) - !! From a set of envelopes check if they are the same line - class(injelope) :: envels(:) - type(injelope), allocatable :: clean_envels(:) - - if (size(envels) /= 1) then - else - clean_envels = envels - end if - end function - - function same_line(env1, env2) - !! - class(injelope), intent(in) :: env1, env2 - logical :: same_line - end function - ! =========================================================================== ! Initialization procedures ! --------------------------------------------------------------------------- - function px_three_phase_from_inter(& - inter, px_1, px_2, del_S0, beta0 & - ) result(envels) + function dsp_line_from_dsp(& + inter, pt_1, pt_2, del_S0, alpha0 & + ) result(envels) use legacy_ar_models, only: nc use stdlib_optval, only: optval use linalg, only: point, interpol type(point), intent(in) :: inter - type(injelope), intent(in) :: px_1, px_2 + type(envelope), intent(in) :: pt_1, pt_2 type(injelope) :: envels(2) real(pr), optional :: del_S0 - real(pr), optional :: beta0 + real(pr), optional :: alpha0 integer :: i, j - real(pr) :: lnKx(nc), lnKy(nc), alpha, p, beta, X(2*nc+3) - real(pr) :: phase_y(nc), phase_x(nc) + real(pr) :: lnKx(nc), lnKy(nc), t, p, X(2*nc+3), alpha real(pr) :: del_S real(pr) :: z(nc), dzda(nc) integer :: ns - del_S = optval(del_S0, -0.05_pr) - beta = optval(beta0, 0.99_pr) + del_S = optval(del_S0, 0.01_pr) + alpha = optval(alpha0, 0._pr) i = inter%i j = inter%j - alpha = inter%x + t = inter%x p = inter%y lnKx = interpol( & - px_1%alpha(i), px_1%alpha(i + 1), & - px_1%logk(i, :), px_1%logk(i + 1, :), & - alpha & - ) + pt_1%t(i), pt_1%t(i + 1), & + pt_1%logk(i, :), pt_1%logk(i + 1, :), & + t & + ) lnKy = interpol( & - px_2%alpha(j), px_2%alpha(j + 1), & - px_2%logk(j, :), px_2%logk(j + 1, :), & - alpha & - ) + pt_2%t(j), pt_2%t(j + 1), & + pt_2%logk(j, :), pt_2%logk(j + 1, :), & + t & + ) call get_z(alpha, z, dzda) - ! Bubble line composition - phase_y = exp(lnKy)*z - ! Dew line composition - phase_x = exp(lnKx)*z - - ns = 2*nc + 3 - - ! ================================================================== - ! Line with incipient phase gas - ! ------------------------------------------------------------------ - print *, "Three Phase: Gas" - lnKx = log(phase_x/phase_y) - lnKy = log(z/phase_y) - X = [lnKx, lnKy, log(p), alpha, beta] - call injection_envelope_three_phase(X, ns, del_S, envels(1)) - ! ================================================================== + ns = 2*nc + 2 + X = [lnKx, lnKy, log(p), alpha, log(t)] + call dsp_line(X, ns, del_S, envels(1)) + + X = [lnKx, lnKy, log(p), alpha, log(t)] + call dsp_line(X, ns, -del_S, envels(2)) ! ================================================================== - ! Line with incipient phase liquid - ! ------------------------------------------------------------------ - print *, "Three Phase: Liquid" - lnKx = log(phase_y/phase_x) - lnKy = log(z/phase_x) - X = [lnKx, lnKy, log(p), alpha, beta] - call injection_envelope_three_phase(X, ns, del_S, envels(2)) end function ! =========================================================================== end module From e205548560ea0866b95723a8c30d2fac2b9a251c Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:20:38 -0300 Subject: [PATCH 03/29] testing of functions with numerical diff --- test/test_f_dsp.f90 | 103 +++++++++++++++++++++++++++++++++ test/test_f_two_phases_pty.f90 | 88 ++++++++++++++++++++++++++++ 2 files changed, 191 insertions(+) create mode 100644 test/test_f_dsp.f90 create mode 100644 test/test_f_two_phases_pty.f90 diff --git a/test/test_f_dsp.f90 b/test/test_f_dsp.f90 new file mode 100644 index 0000000..b55323c --- /dev/null +++ b/test/test_f_dsp.f90 @@ -0,0 +1,103 @@ +module var + use constants, only: pr + real(pr) :: x(16*2+3) = [-1.9505250424097711, & + -1.1667253258630770, & + -0.50083650782086886, & + -1.4403476424022568, & + -1.0811890481373283, & + -0.93383702132049806, & + -0.90796780830843316, & + -0.78367628417639512, & + -0.74775960561138821, & + -0.71024720094329619, & + -0.69526113124061351, & + 4.0543875795062009E-002, & + 2.7816909475822125, & + 3.9157194615109057, & + 4.5068354677502320, & + 11.360614772439938, & + -0.60278167802332649, & + -0.33378919660713791, & + -0.10052077666369962, & + -0.41963150416019224, & + -0.22595193414305051, & + -9.6633793016200953E-002, & + -9.6893816488812616E-003, & + 3.4503360665734797E-002, & + 0.12298863627470356, & + 0.14608613593991784, & + 0.25088155199109413, & + 0.92977835950126353, & + 2.7781340724066719, & + 3.3719143186445981, & + 3.8606337057567375, & + 5.9948813972520849, & + log(232.95614892378970), & + 0.0000000000000000, & + 6.3025709559471332 ] + + +end module + +program main + use var, only: x + use constants, only: pr + use inj_envelopes, only: from_nml, z_0 + use dsp_lines, only: dsp_line_F + use envelopes, only: k_wilson + use io_nml, only: read_system + use fordiff, only: derivative + implicit none + integer, parameter :: nvars = 16*2 + 3 + character(len=*), parameter :: infile="test/test_f3_pt.nml" + + real(pr) ::F(nvars),df(nvars, nvars) + real(pr) :: jac(nvars, nvars), numjac(nvars, nvars), S, diff(nvars, nvars) + integer :: i, j, ns + + call read_system(infile) + call from_nml(infile) + + ns = nvars + S = X(ns) + call dsp_line_F(X, ns, S, F, dF) + jac = df + numjac = deriv(x) + diff = (jac - numjac)/numjac + + do ns=1, nvars + print *, ns + do i=1, nvars + print "(I3, x, 3(E15.5, 2x))", i, jac(i,ns), numjac(i,ns), diff(i, ns) + end do + end do + + print *, maxloc(abs(diff)) + print *, diff(34,16) + + +contains + function fun(x) + real(pr), intent(in) :: x(:) + real(pr), allocatable :: fun(:) + call dsp_line_F(X, ns, S, F, dF) + fun = F + end function + + function deriv(x) + real(pr) :: x(:) + real(pr) :: deriv(size(x), size(x)) + integer :: i + real(pr) :: dx=1e-2 + + real(pr) :: fdx(size(x)), f(size(x)) + + do i=1,size(x) + f = fun(x) + x(i) = x(i) + dx ! *abs(x(i)) + fdx = fun(x) + x(i) = x(i) - dx + deriv(:, i) = (fdx - f)/dx + end do + end function +end program diff --git a/test/test_f_two_phases_pty.f90 b/test/test_f_two_phases_pty.f90 new file mode 100644 index 0000000..dc7fed6 --- /dev/null +++ b/test/test_f_two_phases_pty.f90 @@ -0,0 +1,88 @@ +module pt_2ph + use constants, only: pr + real(pr) :: x(16 + 2) = [& + -1.9505250424097711, & + -1.1667253258630770, & + -0.50083650782086886, & + -1.4403476424022568, & + -1.0811890481373283, & + -0.93383702132049806, & + -0.90796780830843316, & + -0.78367628417639512, & + -0.74775960561138821, & + -0.71024720094329619, & + -0.69526113124061351, & + 4.0543875795062009E-002, & + 2.7816909475822125, & + 3.9157194615109057, & + 4.5068354677502320, & + 11.360614772439938, & + 6.3025709559471332, & + log(232.95614892378970)] +end module + +program main + use legacy_ar_models, only: z + use pt_2ph, only: x + use constants, only: pr + use inj_envelopes, only: from_nml, z_0 + use envelopes, only: F2 + use io_nml, only: read_system + use fordiff, only: derivative + implicit none + integer, parameter :: nvars = 16 + 2 + real(pr) :: y(16) + character(len=*), parameter :: infile="test/test_f3_pt.nml" + character(len=:), allocatable :: phase + + real(pr) ::F(nvars),df(nvars, nvars) + real(pr) :: jac(nvars, nvars), numjac(nvars, nvars), S, diff(nvars, nvars) + integer :: i, j, ns + + call read_system(infile) + call from_nml(infile) + + ns = nvars + S = X(ns) + + phase = "liquid" + call F2(phase, z, y, X, S, ns, F, dF) + jac = df + numjac = deriv(x) + diff = (jac - numjac)/numjac + + do ns=1, nvars + print *, ns + do i=1, nvars + print "(I3, x, 3(E15.5, 2x))", i, jac(i,ns), numjac(i,ns), diff(i, ns) + end do + end do + + print *, maxloc(abs(diff)) + + +contains + function fun(x) + real(pr), intent(in) :: x(:) + real(pr), allocatable :: fun(:) + call F2(phase, z, y, X, S, ns, F, dF) + fun = F + end function + + function deriv(x) + real(pr) :: x(:) + real(pr) :: deriv(size(x), size(x)) + integer :: i + real(pr) :: dx=1e-2 + + real(pr) :: fdx(size(x)), f(size(x)) + + do i=1,size(x) + f = fun(x) + x(i) = x(i) + dx*abs(x(i)) + fdx = fun(x) + x(i) = x(i) - dx + deriv(:, i) = (fdx - f)/dx + end do + end function +end program From 289efbf1a591da368017e753a8673a89cfa20d29 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:30:41 -0300 Subject: [PATCH 04/29] No longer using this example image --- figs/example.png | Bin 38543 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 figs/example.png diff --git a/figs/example.png b/figs/example.png deleted file mode 100644 index e5528fc740292204a32112f3fc6f26f05096ec06..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 38543 zcmd3OXH-<%)@31z5hN)pK|~BFAUP+oR8TOH1e7dUK$3tYQKS$N1SMw`NeW0tf+!Ls zNR}u_mLx%d*#x3>KzAm5C{ZnMFrwD z0%6lU{<*Sk6aLHnnO#2kYrBntrX7JmQ%?Se;>AfCW&(koph!He>J&B6<9x&a*M{u$ z-M@;d!ky0_TPAZ?xAil84S8$%vP{o-S4D})!wI9v)B9d#p5I4IG(BDO&>)GC=fV!N zL$}||{T*0$t^e&h{yfC4@bXAeOu+KXcXrW#-Fx`!Ka7`5jRy_PIBmuyu%8|5puyh+ zc7r5q{3X{JS3-fm%^s{!;jdTJW4IPrpCsrnwkZP zRs{tGdDgx4RaHx?%XS|eCoON^{yRUeQKR47+?=jk8WR$dn2=y#Y%Jn9asBq~+x+EQ z2?Wp1Dg(E<;kro0;G@z`)1IXFSy^3q*1E-)4Gj7_b1c$t6kOBPMt(ojDgXJpA_ZEKK!Cn$)xK-ykn>i^9{&fz0T5Mej$C3O{@p86Fns z-SO7CK1y8BqAjhd>DrMaN49SDiIKi{;`s6MFJF2}+#*R@c~)!FKR-WZmE2Ks>sQn3 zdkYhX5-2FCcCku4uNnH0f9J4_Obls#X<+y6-2sXDBKRJ?X`^bqpOKM~WV*UEQ#)j0 zVxsXoBrGf`Ir-PrR7^|^Z#&z*eaBCox=a)_sb{nd@be1}4nCZKE1zRNdGciVGiC>e z`P~O)NWUKA^Wr2l_m}>haSA@FTxkC*sHkYQK3bBMg(dCn+vMTDDw>*_E-o&uU2iVO zy085CwJ_O*@0(*DsSUTAXy+Fc%zF56bE3O}>+C*`{k=a7_!i$LUWz_(y_@z~N=nLy z50Z2@OI!>1f8vD3=@`8(ImB`5RQl4=($CL63|AkVe_P_ZEKhu^oj05}I5c$2$cVPl zlY%m3@t~9wADy)8()8TOH)_elGBvfvb!KJW zB6EZiZ{-N8o3}ndK37xS1$*}F$<57WII7f_$Fz4Z{j;JXDWA{D>+kmRXz}T~uKY=U zdaSe5W25b@wl|yd@pXDH2i=Q&va+Q*Ys}Moqc*(>Tv1WbGl zir-dUR3U92i!9X2xA}?7&&kfd6fJS2T}Mmn0#VI)?;IBw7r!n7;`kjUOiY-x$Y~yu3vDt+uq$wjf%U7k)IH{_ zgO3qkB_Wh(*D^dIp+8*E#QcO6N#X3-TMcoiiS{G)(PQaU zTemVi%gGVWy0ZtrAw=Oiy5nx&$NOLWm`QQt z-*t(cOH$(a4H9!qEPa1oue)}>%^|wm?HTds1K50heOv$9h_M_#e0X@=x@6h7IZ=84 z{{0NcE=k>;V##XI)za#GxSip`g$tc|*0T2&EHkNc^YYTayUwpHSVyR|gc|91Z`R2- zT8|sMxG7y;B)w4II{DQ0q8;0d#;A(U3T)(mA*wIzPUkpr;%i%5hiB-ei7wgaOo4UNc(KF&(Tbqc4L}z{dJEV!mm*;0@W^VB2 zd2Fo5$H&i~u)+$&%^3YzSKU}&v*0oq%irO!XYby<_#t-O%*+g7l9-q%AtB*kwGHp3 z=6pHtzce~J8kdF*)nl5vjK4gyPrH%htBGpy4dUIWZsaIE ztJ@hG6~MROP3A_BxR}_qQ&B_1mD~sTotA(>*M~pX8S$?Q%q^n$i~aVMUi|exKXjh@ z|MJ7K754_4gEFrBM6FmtYl>Z$ZBazppPnK9S(+J$ief>!JpF*GwdrkszL=C$6sgc* z+}PB#UG*34*H>?q>Hse%CnpC-f7yerr%#`@G+9pX>(fJ;M7rpm1YW<_Ay1rgSRMbSW#^cwoUmYDCjTmXSSJg=q>^w zf_9AA^f}T9N<(IJ|3;XMn*`n2bLX%LJ0@+`SDkr7PaZoa>M+)l=noj6ym`(q{wF@f zV`J?uM})A&UJv&>mX?+_Hq*$jk9O^!{nh-&kLeV$>G1GyN^0ujSgXZYOL9d;g)l2( zePMQ}COJ9z?Af!0g@w4x!{?r4WMt^+NW1(QeRlZ#%HnV2KzHZm0~hDEGaRjrl`Vbs zYOpq3(Bk{Mw?FRh@7ynBdZqFy>$1NT+bIEo*{Pln+qZAuvSo|8xp_~agWcV`znflP ze37?x^X4L_nNJ8sq>@v*B~s3_M%S)gyK^TqKK|z@PUVw*S5-2Qf;~Jos;a6Grxq3# zk>TOWks@_1E!i%YhZ#}OMFa#+@bE~xuP!04&dtrC;)U>Nt2fV35c)O-?mx`T%gf8k zGL$?UE#blI2?Xe3yeE#V)^T5D>xBAT;J9eP# zC)+;EG_0g=L@>I#xxIMtBDLsCbu|}?#ENskLD_f7$(paTqN2LH3+zxfy<__S{AxxP z@mgFru8-3E>rF|>xM(Z$iO-$VfR2sLeWachF?Z_JDHMFR0464;Jfyx)pFS1XjnEtr zQn2Yx$j@*8@k556|Ei9Tdw;p_!omXbKDO`Y&!3cAX^bA34-O748Hzl9th;$|&Maq!$${#@12xU8tykz?^a%j8=K z@6EBkzAJ7AX46a8S5SDhwY33e4j+EGy1IHXl=mho8xl5~sMQSx1rma-G%G#*caifP z38mV&Elu;>xpO!JqSinDV3+3S>%M+1bYF9sPj#B<|B|R2NzcF#9Thb(F##yBd`6Q% znBHT_bp6VULB>*DwVEYry8Zk16}kOY&)E3vvpf92P^hjuouzKJDIFbZCotHDkmlv( z<&O4tZtdKO)Upc~FKXvnDiTksrzOo^L-6@lwYKJ5xQFcgBE-wf3niSP5w&@9S9sYT z4uiIJI@ju&8b9Y>tIKn%6=`YQ#)H*8J;m?eU(47SstFBWM*(>2z5|(0O|2CRFrNB- zCxQJC=g8q*1ETD_sT%2Rjg54cnYp=aIQ>pePOh#}rI~4Iwcox4kV-bz?;Sh#MvG>t zPPEr`u1?fJTv$+QV0`=@>d)=2+!>#1{R0EB()Z4Niir~7=TB}*NKQ60F=3*m4UCH7 z{=+TlIFTekZMr-=gxz05GV^6PhKv+;PFJ@VwVcMsZo2Pu*=be|TL%YLS$9`CxsUfL zx0d{!`f#H16c5k*%*=ji7yFyM1%(LYhK7bm#d!sK>ql-$Zdv^vDQYbsC^$Dei}O1> zH@9=g4xD=OqLN=HPUM=mrkp`2NID7$38`MWq9=JuLA(XGTwh<`%B-xc{Q2`|S@+c{ zYhMycuVl|;vr_Y$%Zu3Gs;d0Ezkl?VJqSNHJp4Q)q%G;{%Q4$Pk=vcLjU(S;Q%j93 zEpND)XQZT@P|-lHnPGnL;DK*dW8*V%tvxDbDABjt_eMr+H(ECfdjQS-#*fJ ze_%G<6}ct3R8BD@HufMJ+dB5Ey}cc0{P}Z6^aRYzkw_(`C{jok_|uGhkdur4y%G}S zuU_Ssk&%&>md3et`@4Vu5fTx}&CLAuJv}xihOLWfa}(7(HOV_6A)&8X(G@6`SiY)W ziS(mak|c07@mW~q%ckb$%ds+bEuvCvzY)bmA^lIA*Hv$wuCA^|j=Fr|g3HFb+tsUA z6&1JJQC=GLF@19>`e95cmkU-GubiA5-{V|dz9g5039f9}jTGHt7q!eG* z$*ZNJmLI%%3r= z;3&tgkKEQ-4&!Y|3~XXHa*B$-eiYwR;`zcIc=))8h~v*s-p`&rTixsI?2Mxi5FYU~ z{oOm>?3)@Iq00#w8G*^xI5{Lz_-<9!QfVCXuNUPj?-3$&baZ|oaws=#l9!hs{ubMd zJmNez+^YKi`!h5*SZ$Bk{Is;R%7BBqC9cksojI%$_BMC!;G|~z`}-g7e~`becsqK& ziq+;yiL~Fxe9a?xzBs;fNiw{1hFE*J4YA1a?W5BB<>{Z4c0N_Ty|J73iVgeu`?Crd zDHs{O^7XAoAs^`P$GR8Xm_p-|s1$y!xeWN}VQy$hh=49YE1#~ZsVObm&+vxz)rGNu z2+EA9{x5zg*+PPXeovklg}mrQQp71^*6daVOcXV5NgAjO04g{r>wZH`Z60kOsv5cF zRNhFZ7&l*u3JSb@@=31(%W%p3BWmV$nkHhHc-0Nw#mTO`1lCaF+|b8+}n~JYkEr|8!oX=#8|rw7gEWsy?T^>gbr;;Al)zIDPuQEMGH93eNGLKlbr! z%}q^9J!4Kc3H>4A;m%9A2A|!MVq)CaGx|`Wl${+>F7V>Gpr9aTJE%|-!qd~;O-@7F z+FglHD0GkfIqG!%df3Aj(S4;9Au}FysBTQlW#NaSP{dK^j|&P0Ka>{#(j8)Z_ip;N zZrVn~p@oF4fk8rpbEUq%KGDd^s{UZWem*?{UF)CY%{?FP`YK$#+Su7y$g8oLT^)GU z%E}7O6hcT*Q4xIzPOtt+rSI=2b!3?Umg8mVecyQNJatBz$+wtVwr%Lu z-{EMyHbK{_s;Y{k)U&((NgJY^MZ%s$D!lt!$8E2mU`N4&inU;PBR=`zH$5I+$bjukE&Fbmn($}c!yinIw1=KoGka! zuu8PFwY4D>GcHfsJ2-^@R5A9A9uxHY&DSt*BqO&p|A7=ieE?P9u?$!PLC&KEmes5{iX_4!sB7zw^R{tt&0?mmb}ZxZ&3b9F&Ozt z+-X`5pc!eT)&9DMh6WPsv17*=7#Q#;^6So>e!vi1$B(l}I;MTC&@L${5fKqVT-e*% z3J43cidbGp2SL4S*ZTTjYalGZwFLmXDIAa!*1aXs0=K>*u*ej|?;q0f7V3%z0EvM1 za0a8LoNe1P^ihH=k8RrSy+XzL>D|o_n>K9%`usLz8!PLvj!#0J#wD#AZ$k#Eow|s# zR?D{2pgMWOzo z_Tai}hiu=zeT#qZ>FG(ew?bTBzut}*#_l3VEid+ZoIih_;aOIeovJFM$7@5+#B75v zJ5HWV1>~zJFBkfd`~mBScg_MrM@L78a%+BaWCXQE83{$vecMq46C+bNJ zOQI&ZIr;_$aMljU%F6l{XB|}`LzsF0@@G+{EbNRby-!#r_vUTHi_cC^JAsBMEiDyN zHao?|HS_0>5-QrZZApLowb^Nc=I9>%RUr@rega9c(3p%i#3PKGRJ&>Ed_H#+*uC-b z@d0%qn6bliehtKA?{3-ExL2=m%^!6&wO8^B3vr{DqQv;w+53QiEh{NP^jE0}auLN& zGYRMiie2Wxw3IS3iMnFfSj25#y?klNGp=!;8v#jR&jGduB4eeepIuw*&5C~gODS4n z9+x7!I{kTl?Qci@rZn`QS7gzncV(NY%)7A>IP0>_ngs*|f`Wpc`1^O%H?^RR-dG-q zMXI40ucO}01z^53(cv;JL}l=N6Hx_=)PWp_$i;r460 z2hf*ZJ0c;Wb>_?fwnIRto3Fd42Q4~UpT9UMS1zN!R;v71#wRZHuuBklHt z^z;oJUH#9GZ~|&$q$4NA8i;`V1ok5$R^3?ho}{GXxepW)P*TpHJ7?9NaetD%fk^Iu z6&@gA+7)!mt}i{9WDgo|&os2PwcYpW86BDLjVMM+rC||k0FK2^t=3lNfp1tuEcbPt zmYd+)pd!RMGmF{u0c^9dvUbokN!3OO*L?i=v9WOm9GayG1;0Gjw&BUI(NQ9iXg}2* zKFMC^8yu{HE}EO0J4VX+ljB}mrta8H1k;cmJ9nZ0$I7@p?oXsZFXRAH5do9Qy19G{ zmL_Nln}K_2Zfw*9Ch5|YyN1;Uh6fS^Z<-PPx?rp$YvlFCP+%D_i`cLPpBZLT68Owr z`vr_p7QhCh1(yYqx5hRA38ho{J$^ja_BPVTXKO&EMTF%Os?|}i2M>yiivw81Iok9Z zo(v8TYhAw%PBkvy>=hEmk0NIks>SrcE$Gt^0f%sKaIi@^(K0fYJBDlwplAWb;_@VWHd_?GjBEaD6+uMOH!&Xz*EP=FH#4-yJj@YA9bNHb3Qe z?%bLC!1H+Z=xQWT6-pi=0-VZ;lPABo_}2_kZQ1gryqs~*9&ATDn1S~87TfkefBuYr zPX}F>oRkz69`4U7sg6Eu*Gz|rfx+x>9W&}mVq%432o=iqK@rO^N^} z8t%%y4gAwi!{NQMFnO=E`z+JOodX!5W{bm1|l3)%kVb z`K@Vu{QT%73?#i<@Uedpah3kecb8{x;aj&=-Ed<#aA0g|YP~f%8XU;er%zFwB2cfU zrf~Hr0Kg6a9G@Ri`#pX7ZnkC-sdaIp;|vk(`aTK!;~;wO44iXyT?H0LyDcpxRqmu{ zSyflJhRon=ip~?z((r(xp&|16)?NGfjH&{G-)|N<@uQ;=71g+TvlvwjIp^&ydR7`*q-n`Kg|+>T2KKWduI>sf^6bJ2o~+?jbinII5|tKH10rsi`T* zjBe%SfS<4L%;MtZOPAEIT;UTIjuyEc;qAR;IB*l4!$kWGvJVWkO#r z|X5MQOHqWFKcMy=tpD>Vb9?-dP93@Al1 z#U8F!eEoXCc<|iKn`?lO4!h0JeL#muA2)XaTVd?(i9d%$SC*D&8_{i{fbM;^1yXDQHAA*yEQ$P5Qbh}TkZg=ybf5q%#3*U;@Vq@S z3G8tY;QR$2`Z)>;3T9C&&IAB}oeKalJL!0C7OIYt@g51)?!O_VEc|lP)YO>aM972NhE2GM^F}8VZ71Q?mm{Q03aS z&6410y1Kf0OWd^41R;RHB%_6Iha81(nE3g#Y*%Ylb@fc8Wx2G*qx}}!ckkA?a%E|z z@*tLsjgc`M@BrviDMo4;wXt7v{*(9409MI!LTk8w{tn90#v`qQYy(S?Zn^b-g?>Saxjgb`b6Z=9!}#8SN}O#2Uft0pYP-Ll zpMNqkF=;se!jjgSOrpzN_D9$ev)$KL#(WI+8jYi!@3FgZX-3;mtHPR81u%t0@mj+% zu3GJe(7;r-aMo`x-Ls&eCj7ucJUTWu*5AMG{jHjk18L6&L~+MPKc+@Tj1{}^|y3BE%{r{{`tcwF}=Nicw&NO-@Yox$HV7AgBh8w`vadw ziQ9!~J=7Qvva0p*N06tan7;@GNo8ndMJ511YJiud2Lq;a?i>{z^*+A1m3XS%Y_6yx zKD!UPudix4#)N7XyW}9d==tU3%*UO6@agks@OI1o=%3x)-HCsYEq2zMCQn6LmzYJg za)m`s0rtMNEkruA{81FX%)sPBuDc{LmU<-j)S_&ZB>V?p(6o}uT{{mad(f8uh?$sZ zr#W;9r5l9Z--Rw~FeMO6!j0;6U&u<}at9-J!dxtX;@Ai6+S_ z+G6aVj&`LnGwZJ!j`uejdSn*DLFl<*r@@(tX+r z(u|Kcw8mP$qL@ThCzONVI)ZYn9axjTw~)4;0*IBFnHjd%-)(vJ2wiJ-QCLBNGiBRg zM&7Y|e=ELhT987-fuaImVj@XW8n6*V1A-D8B58^jf>nyTlcjR~ zzr6r=#_n(U91*K(QE(`dICA9ma`aSaw+(*NxSHsx2WdFHi78W^;mPj>9sV` z_XmP;>!Eu1`O_!PS0!aj{HT2M*1 z6?-Ei&aTsOY*%`>#V02#Zn}FiMXxhl#!GkGP0?r%T6fXDfy;KtNG77&KexAsah>s= zlGagGeI5~^RUBP#ip&1&1o~}>n26=8)1f>`^3(BxfmY~Kv(QkndR%>!$avtuW7*K* zE%W}g&>lGw&>+XB8{fM1>|?opqMh)Hk3+#;wMHQ~N4slmywK|Zc?s7_@BLGD>1YdJ zKTaY{45DUDoMj)bg5185WB#Hck1YF>5CqUaYhrJ!=;)My`P#|E+rmGIgdSd;1SV=0 z9tZb(=-$w3a}FAdkXmuH?*gkpTgiPO=N}an1&CHN5xOa8HUgjtLJy_-qkI%=b{5Ce_BI>EE#|aX#4(MR9JW@ z-M(_jHhilhLchPif6QqRet_guJ}_QsLpPzJn22*|u^MZ?Sml~25 zszmbNy-RRlS?;fYbU75QyyfOZHu6&ftSD>#8Yu3V>1ee3O&MT+5?<2-gu zkom%N-XWCh-D+JPkJFQqMzImwckaBXprE&2FQnw?xPav|G%x^F+}rm+PO$&iuL2Bj|xq;(|Es=tR0xC2srp{g5zB}`4tO%9vpmz z2*_B9k9e`gn&bZ|O;Zm)(qR3M#!Ho^v9rGoWGo7EeXMM3DfKbA!#xFdmEdoey>~=B z^WMo65gyJjdHex;_42mQYg|8pRe5!aMu&!4(Ym5t1Z@j_2{8*N49N4G6?1nGjmQF6 zJ;WofWbC5(j5SN_ZyJ6l_nqZ21ifOk)wh^Hp7SAF?%gpV(9 zlCeQf)T(clkR}U(J8&?viov$+fM}X>-QTd1|?vfcJir+t{hySAM-g*A^KTmchE2D0lYk z5TFABUA5UzdFXPD1Dz!p@L4MfIAB%G}IMxa3{luG4)RXUI$fBMqBW)lZwS zNvG2*WX4roT$8?Satjst%Sk-|NvL~SSy}jG1ABZju;D7##ckwICH25bfn_Rp+(DZT zuM#vJS1gz1RtgTXAE1ADU;quBpx#FUTIhBj&G(%A{7`mcCfPkt-0tYyB-R7>q*)vZP6@n}db-Zrw1c z@-oVFGE_+<;_1qR9;;|0+Uq$!fC(b)46Oll`BmcfH_RGYm4}-c;xt-#95C)>VdXLayuhaQD=MwDv8~^}eKV2V zmXQpq1-S^On+bHOSR}vcRYuRzI|beUS`+ZnV0GsE6HBl#m98_mUm+0|QZM%4a+FvZ z?z{Kvi<2{HV%p@j)v>X$A-54YRm6)IQ)jleG&dWT96xcwe6_LsC;EJ{UhyVy8>9^6 zaI)8=>$Ioa*w`4f<|XLTAXofsJa)+a0W}qV>geSd?jTBRRLHVC7lwxnvD%4NtjgIECs4Bv$?2rDmm z@3?h}!X_C7Ahvjp@_**x;Iab22BL0!a#EQHwlp|9XWgnhzwAl(^ox*?uM_?$&=|H- zQSoZ$K0%Vul#-GHY{=kJV?b&qG-3tV6VM_F^%bKu$*%oj2&g>C%Zu*O0b49M0(xN# zK%axA)6{GH)b}69ck;d6XZ~r|PqZaIlk5VUK|NPv7r%X5oSjW+va+^5b${`wToPF( zz6#oh@(1(5gStbDQV2N!6*F|e{vI?J2XEZC0gQ*Vw9LaoQ`ylG{zbu3+5x~$AS zl-q{-`apBr0Zl=w$%;N}1U(d4JXI+1^=l}@S#)nsw1@_gDq*WF8uL$kc%RDj`#V;; z*Ur9Tgr?&3uYTs<1{N9e_w>D-p44~dI+)Yhiy z{GRMB^_W{}7KNAP`0;JKb|K*e*F9tAh!C}A3lYEH31?3=s7HWNoDqFVKht*53te3r zX1jLo3=UkOn)9O|TDK%!t!0ANY#^!_jE?VFNB>cjpF@X!!v1(ma+^XSH2gFWB9n~o zH<5Q=Q(Ifm;v+d$B_%;Zv+KY`D``8p4b9>HhPF1qm4mrfJ)o?^kvQeJ+tPI}6A^V_ zp719K`#a*GM9>mFNnfFK$h%du7gP$GLPy$PKIHe8tSx9G=Wsms_5IcO(g$V`EX-0m zdOO7OUBIk~Et~23W}MJN4}iEv$*J2ubxoKW`Muv3ZY0etW4f!x#>23X*xTD{<(RiZ z<`5PZ2KIpn(;`~e1CrfiX<8ouMlj2$ItcK2WNC!sn<<|;^7|0dY!FH-!93_^SnHh`P@~pn#+C3RX0Nr$kvFHOcqEgU`9QyV@@u z`e)8_?Sb!#Y&mh9>;yV^)qs>rc5^%K_-9VtfK;buQOICsCI1ZZvR?$VsajhTHuUe|Qu#8V$ z9vHJ?W)Btuq4CqHwgw^fcQ=+l@0L}H5+l?037PYCzI3OdYV8Wjw5`5 zWTYL&7_lBJ^S_$U1HwQF0ZD{?pi1cQ_*zwE3UQtKo;AF4xV*FUlmxkKbwm+}Hzevb zJ*u5odr1#X$0E`1u*$gR>;ERNFRngN?(S4KAiDMRDQm*Kuw^^hCw%O3%=PQnA%1jx zUNs;Pe)7#xMZ^>rufeQ?Wuztwx*vWEcFw=pZ_mfgV4XTmoKN6oxcu*08IBb zE$u)!F&pV$BklSf;N#g@+cHNxySbO=1K=%7y0EddVA6!2O=1(GQQT^br1 zaK9kC(5xek08rD>(Sc|ix_jyJ<-FY7`kI>NmKOWFcSVJSpf`dGKgG)n>H40%eKr_~ zp39SceIEyAkju}yUxm*b?wG&8^kAv$9ccs40y{wibNh##4#<1Sd~)ItHg7r3F}tv^ zcTzvH9V}<~!Fn`qTe^CSoL{D=UnWBS@tQk=+9dmZT#Ha983aCUc%`=x}Ao=VS|mp>FTc+ z$RLG#eFhJzTs~ii66$+gO6o?o&us&+?4*ZQjNR{VU3A}C(d)5MH!uTI5;Q`JeGI?h zqoN{Lef^&RS>@qO1fm(3Z-7o?$t_hvOiZCD=1`9^SzqmExZ?(E0XOzLFp$=G_JUw5 z*l>;uN=iyNB(v8L?R48^_RYgC$pOaWnwlDPA6PHp7$q|9-Fvj{BZ2({cAAZi?X6ah zf#kL-iCsU1f=J`c5c=cd$URQ}hYttEY@qrQvn#C+z$rNaV2n%*q6qk^)C~X%+@adF zYn_;nD4oF%i7f?oCh}dpB?0niZX}pJ8T26zmSd zGzKJ5Ufv?;aag)g7_hKB_H|)kTB>R;_HTLsTL5GRyb<+~ybiJ+CBebXNLr8r0hl*d z=PsA|y!l6h;ta~~{?h%)1IVT4TMS6-*x27gLu!@lwMOp>4<00#zegQhg`Ch zj?9F?i6fkx_LH42dN|=O#mnK>ew={l`8sG0VPEtjKHyTO+B45wx-T7kFwyNbDQ+J|nud2HG zIs^jPqRe_riotxj)aBoJSM}|O(`vYo@TBrLrTxP`2O)6M!KIAeKk>ph@`%pWX=htP z;BhH7YvinbMHIfbtgU0i!>gr!Z*n-cvSJKX1W7r7L;{6HR*yywC``=8Aq`#8&^VTW zOrY)|EGhXHDezCJ@KayY|BZEPaz_78Y(aoPY4j5nA_y^3STOXR`NVB zkPoc|cod*}r7*txX0U@?FLn{8n}qEOO_@>7aq+cScI47{2xPjf!JmI zP9Y8@LcouS;TcK(5N0>+V$c%oH0lHXCjf==o*pmmU-Yx z0xJcxL~R~O1wa699-cL5eP{2P)0ua?AA+L<*qY3$1qIp8^q++Rs=fXC*|T5B6K(FR z*E29-B`#h7+~eZxoP*&Tlo8a;W5)nj>>%Odn_krFpCABBgYZ0V6{m?q`UvjWF0Kd$ zy7b?bE&N`Q{LudC8o{{%*M{vqbB1h^IZq@SLF?@pm$OYuNfEd^1f&ri9!PLGXA?d0>q=!kjtO!zg9RIb))%Q5 zBI5V#SuxMea3F~}&kll-_=s>jUB37Re2pRzn$bueU+f~V<=%;kG?&Gp8m7h=RzloU ztRgb@7snl*E_!3kxhg`;9Pg|l-=8L;E9&1+H?p=yjfVHjeZ8C9jGQA9X~AwkGyO(G zxWa@(fd*SY#&xqPFNP&x8buL;c~c4}{};|DBtr>@v2xvscfg!LuRt%T8LUF4flxbf z{$T9_Y5+z1pN}gAp$M*koHwk)Zr!-c;Vu(JO8MpvQF%XjU(|+W9&iQ#iDfAz@{0AJ z9SGNw0ccTVU~V&(+$L9y@MB9zu!$2t?-~_t~FKI-;YS*8JGQOeDZ46C>mClP4jB7#2$s^yors zdZ2uGh&y-=lRp0x6Bh@=uRObGt{0kLVGZvM+)QE^S8{N0kd`h+N`Wx|hCw)86=uti z69NrDZJ@27?&2aGsDVa~A`Vu?-oV-YYflgOv~ghvIMu}*C-yxeS#XF4f1`Ny1d$4_ zq&N#r*J;9m*J)`=C?e#F$S$u?e?Dx~dY^KvOf;7j>|N)dvS4sQ10*rP3&b{X#^4*7 z_Ur+KQ9ggZvcCS(pvhJO?bpi65ubzZa(5`N$)S($=vYQBfW&!A(r=LdjvE6V-5nTJ zI66V+qCrWuj*boq5xLV}j!5Fwd_RaZnjkGnkW(%yF3y}jsX%@A?jkw!lq7qTMtW}` z6_Mc6Cv$Pdh+;IKm~umtNbWfgA3l@)?+U(uzx(Oo_FAfsglhy9o01FV^U5wR5_JB| zqVGavR#sN97KanEvxUfqT<&-+Hx<{*xHvx4MNdz#p9;Km-$iSSZ^K{F`N7Q8bXt;< zKx%~n;GF*fa>tB;DiZx88*^ZiM2ugFY40?W41k+9MU__BQTP^~z7VD15 zAnd{o1Ww}RSa@(JMn*Jh3`al^!8UXNwRXyz(lfsYwyMjQzkT}THOWp$FoD%j>7Mr( zV937x`-j53%_}5f~6`~Jw36l-+lpo<}oq#1TS((Q@1sPg(s|5|8@3 z$gE=jy+k)K{Q-JWO-&7x1#q(3gY?CuWH_c`Q&S;>qQRGf!c1=ALtyYlrfQ!c5OScg z!ZyV$V6+VfqugMkJrh&C5n^`|AQg6La$b5l775@BNh~>pCLNvH{FvJO0uVCBFS3e? zSR}7ecrNKpQ@y%=;lg9M1MuFPH*bow)T2QJquEY#ay+$_ za7PhX`ND-yn@ri?LU@BzLrTJ7gGC$pl1!2K1hM<4SYsjb89$H} z6Eij}42Z=1t6zM&&d?u&8&g)c7vk*&G#t{76PuD=BEg~TSA6~quOGqS(6M70mFH$a zsHT<#7;9-oK%!8oNK!QXO^edP`B4_X|S~H&d*%% zZlQ<}hA{2=LNE88aab0{lVPE25-o5ag?e(Xhj{GqOR&e6lW7T;h#-Lf;g z8a)m`Yd@q06nt9DG=qr#`7=Hx<*Oq{wHFBcR~OK+poC1bZ?D>BEYF|S&)kH{0s$ z&`unL2-QdtBU6QgGPG9LW`uJ~uBxaQF6^@Agrm9Q;AVmm3$Tr`|Gv|CR&+531}!Yz z3^8987Ir-7IDz1Ak3$5NDOZJ%p*L~;@e`uk+=jl-W`vm{|w?(h6LLGc%Eu^SVzj_soC&>~_cdFXY3- zNT83`l`RDQN*hE${wPCqM}V$CI%wlOiK9J{Y_f^4D%KyOC|B#nNMxFB=?_REAaJ_7 z$>XA{p>Q1=7#hmDF0m)7PvdSC82$?BD>3Y06L132z=0429ff0MiHH7U=Z+7c-b@K6;*TgYohQ3K25$lP9~#>D1^vCbV%wOx)K$;+v}t%r7iX(G7c|%#N`7UwY89E6p8SDTKxO9tyIAD& zO5|YBP}Pn*2#P{?IM(KNcO7a!iY#1I6t%8V($Wv0Lu`3>Q}W-t%V8xW+tB`*CpXK{ z&;L4;u9L~?bsL!&DTqWiJA}w!#0b<00F#ml!E*(`SCV`oO__&_UCvtfzK?#OhPOkK zIYR^scR^WsPiK5~wmpal;$VM&0s^x$i_()@3^)g9CPg-KJHhp4;T7WMu0z+4^skUj z^XTp%G=ro?#69)%!|Y@2MY12mW8D>=1Dp~N{`oHR=8Z822u7QY_fopP2xb5tXZ2D` z(Cz&V#R6&=+Oj!7kUxKn4gdgzCyNs~Po4VxPqI5atKz4K^9BnW?*?>5CIa;C>i>D=l zI~4<=t)a09=)dao3Ck^l0|#dYIWcr~ec+O(RQSN7oy*5uMEI!D(Sf*luiZ;>a-Tt^ ze99($;QQ(^Vds1IjvYAyOEN42Ejrys#>SYk@g_kiII&cI57K+suk-)Oz(;-*CuGhd3Ym$L+SExWV)j} zxPCfJ52^K9rvmcTy1vZF;JY)pIgk~VnLM1M{%weNp|))8HaJ>HN$@Pn>)!guQPDHX8I#iw1q;h@k&sDR^2fjt z$im3j^7U|~oLOUe(pG~pe9Vg;8F#4i)1Mj9XTOc$nRWdsFPGVFA8*ywVs zdv!+R!+btxsz}HkhE!lm(AK`7s4uzsB5`%K1u^sT4j1P8DJI30{>hP`;<01N{lh-DxPK+#w<-=svh z^{@p*MeEuGhG_A^b~sFs2N9JHUZYGp3|Slhr+Rf{lU$IdW@ASdhu9?$*%+O>=>kPf zO;g^R2fG13M=6Hx*YohRv?1WCfiRDh!C;DqoCHq-wp5p&w)|y zwj&(FB`*>ZreisqeM*b;S!bL;a%wEr&r68OWo23QNEoZ`A&kAgclR+VE`h^*&^@4T zMn$oF-XAl5&@jf^H3j^bIZp8M9_7H3DW<2Vfw^Y7@a&ko{pH)?Q^F|(@aEgC#NYh5 z=IhsoY>x}hXrxCI(XioB9GIzWAe|%xvY|J(zjFs3c^)!)c&-2WxJ1PITiG1^LZysIy&8z01 zB-y-`a6k{o4#*M`hn$_T(YCd})`sJe!wT}}ue7(K3~;1if&}gs63~&)a>y_BTTT4B z>gs#ftV8$-E#Kksqo$@t*9xB7bNFj_6vyEXnle;k3i+8ifWBp+UJo9?sB@MGca2gFogdg8-66mEYoj1m9S#)2{0`7q8 zn`lX9Kz@H!R_IWq92qfs;~#7DR|WjTj!&w;+Q8ofrd*hq`}xla4R^6S+CYpWBT0d; z5yA6v^76h_S2s2{!%~^8H-av@ps+CH#f!Qy{?F)J+uGDGU$!pu+JbHdi>Hfe18oh8 z23O=y7|KHc!kRud1@-s0(=-keZ|HMtzJDdZk2f?DD$l-sV(zO>9AnLIzQOcysCwT2 zKLi!BIQY?6L@>%QJOr4=wJu$QSlxyUn?Zi20Dy9FQ4tux;PCMMa0dnlzf4XpK(_;9 zllrOOgg8fXy|U^K!y_SZ!e5;ZfI?K}NIH6~(qh)Z38)OAKyJ^BGb49})?kFYs0?YtA`&mYBhsV1`l*ZxHgqIpy?{ zdyBjvfICvmkYzztlgR-g3y=jdurRLQBG8H-IkNK2V|3m^fKN~If(pt711l?rC>GKk z1^M_u40)PfrY6+&;ESNxk?${B!alA=wMCS=6{H*kLvV!7kc}9)Z>Sm&>@7Opzy9|d zoghXL1nA!Y6$rmta(T|4YYIGMIbwPOnb{bS44$7xHU@%G8AuD*t*dYo$*mtmdkVkc z1Uwv?CwRJ{Aq8pxJE4b0papFB*P0HR#)s{i7#@=g1jGXEPF|_HC;rL9BOb2;M=QzB zKNz=Cg)qIuhDQ&4ea?l4NnmFBE*|7j^=&7iP8Ea37!L5H*hnw_`_y90a$~FeX`>To zDE8R&v5DOoAmk+FID9W2&vnSL%Dk}0Ff6w3xekgUIqDOS!uMLCM>0?^9TsOe9qV=o_)5i(i-D_ zaHci`?ZRHX{Q5XDJQ$5iH@D&1a6tA|u+I5A840AS|2*v~O@%PsTlMVE;C69|_x166 z%w^o|z{oyR7#tlv-rJ2~(g+JHP(dOHP;gnmPHRg80b@PF1l6+`6s!IE6{Mve>t@dl z9HpU)ZHI6XYt(tupI`!ah_WEs4_47-UQhc`Oey!`L8}l&c179! zKPLJX!i%CHsmV0Fb_e-kS7?nkZz{vmqX3W}W6=t1i6=H7+Q)h%3939~YHr(RwCVU} zxl`Geir(I7XH5Mf{zhT~lS4pMG}M2hw=|Z7fe|Qeh8=n!#z~_P4%L4DpolZ?0OZ8~ zuD>BPFi|rHt#kiQ9A2-$a~McJO5AGJ_<49H!Oe3fAjPeH)4?TdL$muw6D=S_aj1^h zw^B4PDMwt9BaG8_CJe7Xx+nD=lgb6M?{tdE&zpb+E7olxd^H>9CXP>QDuk6?8b1i8 z;rA&D3a-j;8fVdpA6_>G1%|EzQ#fR?7aYw$R0yLCq9U*k4k0{iM3j#Y<6Pk`i!=na zv!H%JQJ|#ONNn?TH*4CwrtA$!f)>0EBPxKn==F{yfK$M@td4s<%=T*6uaEZk8#6B-gjc(}yF9k#t~s9I4$&;Z$-SJXCvYz$phd4h?M@(^$u1A74uI|yzq z`TWH$UnCx?>$xct8glWIqI-)^$v+t!Rmii-q>0}%_S+i5B+gew4+ZbxX9_e2W#G$W zr}stUDr5aaYT0DaQ!4(-amF1;^NbZ3QpQw&7;)2;m#}wViWGT^RG;>MF!~Z^hJo|* zMQYd}p4U2DhYpo>c4|pWCb6WYC;Pp+maKY0jm?OjO0F2V{qad_!gUyn3msbL)KAE! zD458+ZRzfHwRNile@#~RyH_HRR$#4K#>aHxsc&E;@ij?MNOe()$TmFe!Me-*+0#*Lxp)ehgYR#D-nb^X?Q&#)fj;DOV6=*!;~uYZxw>@B$XPli1)sV6Ed3K zvA%wlK@;zMU6+B7JxF#u!h%FTmUq8=dH%dy;^DA9dFdVjor4z#HR$a?OK{UDm4PPz z-@itl6C}G@qc;(f%v8SC7_!RGWt1ExnBLISiwh4APfa~J&Q6ePO;sUo9;c1)nHQ#z z8?YIo_MEb@SM55yZoj|n`({Gj4jLK}KE6?ysmbpqw!GRN9!~xaDv4PZw!mMQ3d0h* zj|$7}X@!Fr`>l6cnb4zerpU!*q#wV2w=Ft=PdFBXvVnI&ox*!2BubF@^2J*~P}8S@ zPxgE^g_BA9=cKh7v(-h#*GoI7)GD}nXr&=AHvB-*$L zeI7y<=hAhqetwi@XBxEfm-(4OHf9vi>nd{6cj&phmrRTUePLp!E=qh3o&ob-Ny6w; zz!6M0Ab^k{`#c7R=zhm$q^{4Y^tvtfAgdgqSHzM1`{&ooqCmF0UuM8t*!LP~XOb$FxisNboQ<8Bc!$kKt#Nx2Isa~NGR+^*#5Oo?rTBt7fPiHXR1 zVaR$2@)hg9*LGM9O(V^w{(HC0G?iRvc%;gO9pkhI$oUQ9m?Tu#3*iiKTEp&YWK^z9 zMTl{R{10vd?qj4iLMt~nX;FieR8YE~EX0ysH838ii^peV*Z|xTh46fYoyt-YbcE~h zwf4ZqFwRK0!XYV{qT=JQjatRdqxr*KYqWt_8}M?syYu-uIOMta+C`!;VdlI8YYndi z99@rL8Ci-b8gA{~%Siu8;^CuYkF{hlR`BAGtJrDGm`^paa0tL0 z3!F-xmiVzY`)G>$2gMZ%c)T$o-Gz*$xs~F7)OO}^HLmZ!U#8Gjl#C(LEJ{+_l%zsQ zvj&l&46&t@L>drobx*W{m$#0{%PB* zwbt|8&wXFR=lWcqD_GXu&1(0+=)}x_cjWi!RrhY*e24!>eg(dJNt`^C5XuPb&MP$d zw8E%M&FNC^rKk7kxRkI^*S=pjQM?aphNq6^2uM`@r%wWIPAnqzX|KM|iIyH}XNdeE zv%AAs!z^4?b#$_}k#%wwxBTTxGp6mw)w0&G_dwnF(Y^GCsC+o?4L-ux@88{rcNZC) zUNo)sr*(UKM1xgRO9o&A6?R~N&iSAJ=^EVd4~Ya=1$H1>)8J=s?b!Sr8kHqwvtLk1ZB;_S5GbI~-2!NjSJ>Dg1Z*YDeI))%_Q zXyoeFts7~~sY<`)>B-59=^^5_1<9^SOG=}Zv2*aOsBB_&dD7Iv1lf}d7kuFb(?W56 z*kHq39EjvH|K`099{lHnQGcV;nh1U{G|_$gEa{!uRFxD?HNe!wL?}jqdX+7Jn5yNH zZDA#z}w~-hUK1Mdwe}sXa$e9>{`W<#wJXj(nQyw1JThw&cb^l8K&g*oJbmMWVkJP>7)5?!-oE^zv3*whRL@2LnQw5Ic5(%CBW0lTf?o z98i`G-`lJYI3&=W&C*yTS~N|6|~IT*D#-BmU}`g#xdTpPBMni5Ux0vRHdy*f zDw*mJ0m;MD*}oc>y!Cei!mUyqC0am~1OFeGr=M_5Yc1mTi`UDx9liE_h!|^pg9LH@ z$bQJP+%7Cki9SsFCkXdsdXLfXzzSr2b|#hxHA0nD=dPchU$NrgtKFq5MRu#)Q4D0Q zhQvnGG5RWTd-1SV`2k9@rLgM>h#*^0A2UyKadez*{8*gP5wEMGqnK0k@YZaNn1fqh z+TJ<6U(#w9mW_bq=p~8tl*$6F;-+@->l0BgTI1nL7Z(@rT##@TOa;yXXUg+_O6cWE z*LbR$GQ);dU>12KTfBAroX(b+?+Xvhotgh_w4o$rHtI~&`CWg+idI(fbOPJIJKz}& z<gx+?!#W@>s90|EhxiosL*`qJG8IXA({0Xc*E#y*Ypu`^8s zS8i|CjL}(J5mh|q{MJR2DKf(UrgDiJ2$sa0n`QsKRrOJMRmElB)dYvaG8(Kb7jBT) ztATPg`w9E75;U5ibWFjZ(=D`LQxrrD^mA?PiNDnmQ4J$Qx;@-w$kdJNNGU13BK-U2 z?OV5)&j+KTj5_nm1_$eN*}gSp1GuT=*GjSvBHn1mrR?Z48NiCcQ z#mRW!gM^P{IYrR;~kVx$^DcNouNRKesKcqaIXb+ly<9oou$ zQdhS3p5;6Ngq5fwCLW&S;ppf`r<^7i8zR8$%;vGAZ8n>Ys+#vE{=&(USE=r0_<1224hlyRN#z|JgUsN3pW z95MSJy$dTIej2K~axq#~S1w?{jH3@7axBfO(+M%Wx^RJzi&OvidK$62=rv~TXKSnM zeM<54rpu++%?7JsJH(* zMIAjo6~*xdos!`+@epskK`ZLb8-Yyo(%gO3sut7);RMbb+3#=ZUda-02eHcUM!fGX zvR}A%Tx+}z->1>{!Mb~&#*P~|Kj#8Z&f|IR6D;r1V-{UQUcq~Cyx?I^g^hXy-hq|B z)%bH*Fvm10kTwP_f|=wI^XF&hXZL|~=Ug^=lfXx&`?KaNF;=3l2u(-io0k*Sx05RZ zPA+N)G!(t+a`uqN(|_-!T}o894Q}hnwy`}#Vzsk}lw>i>8Cg3mk0d0JGP|68`npgq zSfiR<31m45x(9|VS&u09P8{)^jr%WoP*r?S#^~OnIyK^=S4q)wiB*bwn`}i^aXf$h zmF0_RCMWa=o@{>xfnnYz+bSRk?ms3$hQbYm*S!!8fK{wUElc5JrACx2XBg8`P#u*$ftLbaX$C zy*#Xbjo4uyk>Yjwd#uhczSfmsR4`|@+l|BvVtHXE%>bRF=|LPfZY}VmaFfT~RvgHi zMs&8LLqgP1j6?tNx#}KQu(Et-U$0u2@%1b!UA&zFynK zwU4N?ib@@-7G=WQ=5NsJOmf}B{d@`1Bh921_u@+D*`CSG-s{7g`vzU*;zlaM7=e!E zh7N)7fSgvXa$7T3LV`qO@Ln4BWbPA(=5hnV)xW7MP&Ehr(Cy3Pw=y=qf!C%7LI~7X zy4cB0Tt!g(j%yUPzKi!neEpVi71#f{Bs`zrFO>UHd!0ophno>r{I89jiJkg8`)6uN zh>*O*VVs-SMdU22V3`zj_xD@wPcjk`V;FMkd#I@8=H!eI?>_Y+I#=X9P?m%=#Jn;u zRYK~~7lUipt`P{@44ylCc6nryYMrM~q@6~LSz6h)xrL&Aa0mziG*SE!9?^)SR?})v z*ag}}-VJa*T5x|7^*o?1l3oIwpF;bEqlcv7nHx84Qcc}F#5*pbyu**f;6BpQPb9U9 zQ~rw!5ZNK|bWpdt46HuQS)YkofNMnORh2_K$4pqqt}(R47=40xNlXwtzvzZl(>r!@ z4_3EsG#(nR6cx?CvDD536fE!F{jb*&z=Z{+4`OOzcU*LL+n?CBNZ6g5R?q-Rbi?O~ zx~Z!_${Y^T5z--HSm1Kmr?kBM%Pp+}H{lvS-I0F~vEt#%Zg(xa)WxExO~iMW1@1F7 zV0MK3Iq9QgIR(kH=uymbE0vwH#c{KQD0-^D(eTN{Kfhno;>`Hdv_r3eaD{8yn)<%K zr!!r6o7TW_E*rnqS+g;M4~O-dHaTz`m!L)85%_rq5E@r4lK=Yk(aAxN#hlNJui79u z!&33zvb$hP!cmmYzMtGhM7cF-Gk0G<9cHKNk1PKVVpZM0eSQ+8aKP z@ZZRr7S?%+Hhv@rp*sMdJ(Vqg)zM7}{@>*78>+cLPm4p^)f6v;wcQ*U`R)5tp2r~J z;~rAn6MBc;+04B}+6Uvz!owe>%X}R6BD6yy7*VvXf(9Vt2H29pc?&L;JQpCO>n-Y{ z5kiYPs2&llD52j_vEQIU{e@4w$5}<&XSYdVm*D;imhHcIBEG@}*05QC7tr&R5zb`) zhN=^neKP{L0Dek?(QR}`512ela92J>?Y|!7V>L8i{ge^~)dRATgcE}Q-=0*9Ubr6hsJ|m#wz#_ zivfVFTUWl(MU+C=k3MlCHn@D$rM1@adZ%V;Z|yJ*g?&$ezcUP?5$^vHkUO7S0Sf|^ z>9?f^IEVG~5$r z3`#q8>`&ys`j{z2;i(~g@0HhYz=pq)AnB6VT{$wwjQYAp)wh24S)evLxfeJ)3-xkN ziF<`F^o%r{_qbz?h285yft!G8tAcWF)K+?*WBLi7E^QxukD2boQ*dAB2sDjeS2 z+?XzsM4gB^TiBWPk!O_(nnXsbnAm zL-g4O+nx-;PdBB0S$83CfspHHWk9Dugl7Ifd&hKXw{XG_oL-a{7Y0*lO3O*uQQl!C zM7%X8v%d=4d9|xe#`*fETi!@|wd=yDp`1{Xq|rn4>d2H?nVEKlpI;8J9xy$Ii=0E* ziz+o0PmT@3N}F0Nggy_Y&Oe9w6G0i<(aic|3)T1q3@WMbVMvWE7y2#?Vsqi?4O*CE zRJG;AU46Qu_|dv``*%KSZwhCQRB?nsfFtdx^t=p~m%rs%Gj7Z{QrQ;ghRlVOY{$Os zK22uG5J4g~EUee>kDIo|>~p0xc6a^XsizNHQQGZ{sA|u-C zO~o+Lg>`;@^uV z-|RVaj&?3`?_=dpapF(9MfbYo_Y_?|o*X=4YVwc?jO9SqWW%9>Zl=S73!w6R z&Z76~=t)zi^uJHyX1DgYV2kgAdqoRoyd+fp@bnC{vvtCq%o;vT{P|Mot>gl03nkeI zWHD*Fp;flGq@+eh&Zra&O@C{)riX{P@c5xagV+8I?_sFLAFU%bG(OO=AkdcuyP2!E zCYdtE&T^?bPrH7*i#Ve+SI;GrLwN8p+sQ^H_wr@$kCq^{mP`6KTf>@v3a%%0vq!r& z-DOjI$&1XqKFEEqpCdjv=)tV)e8)&fK`yZoOoJ^6g`l$^7jmlcywJu%6qtA{-yBjG6mRR%UU4407__U9LxF}sE&*LU$* z@zuUuhCV6u*I#70|MC0&?b}kD9M*GfhC?SO5$kJ|th~yB0o}8DV~6 z${j6jWS5^%)@pDBRi7iRaD0ki-K8q_LE1iTm!Q;noQ5>lJ$9D%{*wZ+my_4cre+li4bE9>g|U^+>pfBPUc0qHLPm^XuX zcYZi|CjZ~_pDdZ}6FztD%D1zhO!AQ@LeqZt{{1R<_lYnoI0ea)`^{DrN$v_?k}$Kj z8sBaX{R%jSSVD@2Yc)}Gp_#ni-_q6t1jj$2fUa!tZK~}8upzF#Jw!)4`Aq)Ogj5Z~ zpl(-hvSZuY+sVXAh5ME2-TUpEZN9p=#89u*rXvUu}7Hkj#yuTo| zG=KdH9;+gJp*!N3can%@v+i2^@}0r+=?I#>>07TmQljXWX9}WnCPfS_nmj9JE6e|) z35QM?fjbp~Hv5||bA9Ru!mq93GA zRuSzUVG+vug|Do4lyuFqrr1l}s!o;L$AKJfv@n3)4QQWIoF+FYfG3CVxf0R+>sKnq@1)BtRNuIj= zHU?8{x_+i|fSuQyQ>ijCL+WNr?|=VqC2H-iqh=-d?kP}&#H?(pzCB44U3_eg?q+AZcwRcR3lb=W~nk+%fC($kVnC9>McgZtrdS|^x^Ke|&-oh!_`Gb?QP^phf zNQihRAu6M3tE70RtZaG3%I>Z5!Mg~G5}|BC@0mJnn$rq4pxYD8)}N6b9r3=>9Xz;T zNukpM_>bu&v~CH1qA^-&J69MXpW)~E#Cjtq#L#4q!t{J9eU2Nuh{}roi7&kV{>fMf z6~jY6b1stuz`(4iakH1#KFcRA#3 z@)g#i?9ad>cau)c+?l>?))VV*62f1VV&b50tX;eqxO$?lzR1e-{(W6+VCai9~a zO|==T9kR})26fC>fBpBDtpr|zxsS#C^!%LaM@>#|MIPrl>h#!`*)!(PKlQL~JkLmj z^;=F@PQ(GEtn&NsrfxCIhZ8$*$_A@TYx%or{9Wg0#+yLwWaI-4;&8%!hl$d=E?rtz zKI(E{LSSonK!D?kE{aA$-A1DF^_7&0L1lup7+d$$8^Nai9N*IV;L|_n{U;#o@lkhV zPcqnkGN&ESrcl%;N#3(}uj%@3q9mV=sqpB}Ex%uzWb30VC4X~rzc#IJ`pKVZivi60 zYt0&lCY-Q{zd2c2RHqvm_T!UKqJB2s`O2XY(~MT6rmhP#tL+tPlwk)T(T((+5_Np> zJl;(Y$lL1buK>0ohRGM&%N_trNk zsNDx_?U49cPk3DbVULuS@GzhMVpad1*IjfOe#ZNHCrQba&CThyszWF2(9tQmy~kHG zy*s4{OXB4`QXS@efXE|>`hcF{C*dN)+U^n$-W0#Bowa=V+t)_xYabk*aB4zGTyjB| zMQDB9WCOOhx99r2*Fi@F9$S|1P0Bh^Ml;`6Ow!O$3<);fd+hqojF~;3oHcw|UX^Y^ z2Eb_~Lcg)pi$xT=+K3j_tt{csVc5Ept zfv`DrT3mi-qC}s?34L^Xw?}MhLjss}x^wsL!3f%aR}$Kj(40fc?|Ac~$`hJ)uQK{L zKQr-118_3!2h=-ZFg$s(p=AFkj+y$FBDLXdyTzj6r;K`@%)DxDb>ZD2n?*6FUcZ^O zYtOMkT52PMjf_ePx6<{s+0Xer;{sG0b{bP*JH6oNhShJCCwFa+80tbY2)B==x?r`4|?x{F2M zKbfC&ux)x5^U|cpJBQ%Pz}$(jSKR3(MTyH*BL+CH6U2FWXWpw0>H`3yg8DFu4QLxb z+i~i_e?+VY4~ABS*ymWp_?aq*JH!5iV;#)P=^`3Ya_@w`r<%%+%KWrkKXsibkEh4e z)9>r)y}x@`cayGf{cIfcyxx;Vd@=q`Z{ozq077`Z<=%e%Y2JEIb)bdV>1ewm9}Z8@19!#-q2n|{NyRjMY1NsF#n*}1j7lgbOU_*_Zw{;h4&kRkUU zpR}_w96s6V<>f~g#z+k@N=$rod7{7nCn|504p$@$_370bFU$g(EXXo0Xz_T>PFe_g zKuvA7n%u#_g}wS&`}ln0git1UMruS7<@4duD$pjt=!trI%ZE#e9%)U9fACO!Y=pOW zUC|$AbB}%V*%hH|ZZ6sDi>|KREK{36gA1OP7cc61$X!@7C)dyRwDYO;GY5op9M7Uo zB9N-={cYj0gwIxz;OY9XYbHv*S+Jkd^G7Qzt_- z`m|f5%~pycBT@9o{Lobo#c&YHqgI0l4cgi9?tlm2$mo+vR8`f{MwdM6+-SvaV@@6!F|Fi5)S5LVH+Ot^G)i6`(V_A77d)1>X$PC>Tvau5 zC%sv(nh1Mq^2+j&HNAnrP(?F8sHy&|zZ2zu!unf*M|(dtl8K~5HpDe*|9-2?mH``k zZwOiVC2-K>!E;)aO-yd~RoI}CGIXLrz)!1fVPgCK0|v%9JHG7|`sDG++M3#UL_#$+ zC3o%G9~@sh8q25Z!fe{6x$Y? zj_fyNsF{>hQP+@c=ZTSqBM0`n5tnZCsp;(Q#s;GwGn(emAp`h8A~5#XINk2Wo*K%q zseh#$i%HM>g!2 z?54=Nt-W1t<*`z2ylz{8sS-*_zlPZSSZH_VZl1Bpm{grrF+X<~jx%V9(#dCN--e=H z1Kfjx?#u|XRVT%-0MYzUaWsVlB&(oS@S9A7Y+1wfx=pwC`v!|GHq6{ezebrbAkElq z$&#aabA{LXvPT@9oxNR4=Gm#gy6OIq|7rL~d6nHRKKJzLu`=A($s_yY;>?UG^*cxO z5Q&njc=_PDa}#8UU+4PUB8d*gAv~Qb=8L_;>&3i2pW93%@Afy_F7X@C0QcSePKu`lSG%P4QB$Mfl-jH6J85jqbnl-)|BQd7O!2~K-DX0 z+GghIWqUlWH*pnvvYW%ZUSdU)q5tD1)y~}4R3$saSsKH1OqliRx>O?wB zh14hZdpu_CoHN7xpSK{kK?kTx}>r~NDw6F{q&^GL{ckV;E#o0{%C1EC%I?9rJel4TSN zc!)3t2|cIKRrZ&^&__;gR(IebY!5OQVadLLP3vrX06X#_jT8WAX zaD(%y-R^A2{hY1OFaQ?f6W!sg}n>XZ$X%S2MvWy40BoPuhL$LyxjDZoS% z?M#mDtF?H11Z5Y_j7(%m2z32vLJbDK0VF9#{Ae2(EDBWVE7`YCAG?;Z2dfPdp{kDF z+`F&%AR{z-#kwKxK}on09E}{Bs|})(4~; zso1k;&z7q1LTPaLsRf9uw!I~YZTRwvywYtORTNlD#&#rE@X_wbee2Ecf;fbOJkVL- zq!>x68#hK6Z_?P*Vtv6*V;Yith*>gu=9G3I5mn|oMj5tvRmzR~j6=t>pm`?HmY!FA zIsSpG(dfMAB3o%UHpUFRmB?_0F{Z}h=SJI`CkpMO5-uar87b9x`~zu!=HF(%-bjmeIhYs-_n2%Y%z^h!>M(`+(~6|E-;Y#c@M z(*qG-POUr4gO#~Ylm?S^bxp?(OM7(r`rij4XfHh+AjZ20GYG;6R-ixGwa!f(kSXgI zmqpntOUA%H4>b=(@xt7TjVi2*B5QnbjqwTL7NjVeP-&!wOq_&TL)tfepc{Fa!_0E% zRJe*x%QU}U-)I_*lJd?gBxjDU@}S$~w}6uG{1j=|fUiAdIKq~CDShn$1|dSr3Odef zw#SAtDqve~$u~@U#l0<|!ZojVPzq=i?TB`@r@?-m!y)xHqcgHek`C`l7%I-T5UZC; zwAfuhgHDDx-_-vh`pMzx5J$3P5;TT2?_9kNfZfgQtsqk+P=ItM2;)@|*;~1M`SHI? zpFFt;UCUT**s$9XbPQqzWh3k3-2K_fDgyM6JW9qtvfK2oD zV+8k#Hu;T5u16FX6_r|!TVr?^nzyIeNh*k9UHjGFMYW)r7Y`AR#B&UU+xC_>0gIzUPlQ}CX z^Z4`vhno2uTE(bop4;MG9t5n-Y71DnjuTy9zlL@Zvc6IYe-~e7vvK06&jL%kq5s0( zxYb}hvT^`Cp09$*WA+fcEp6@_i^ncO%N8&GoAyG|3Nl&0k5B%he!vw8bBzY5G#&S7 z{akDVdyWDF%)Dn5Y#;_xF@&i>1x`g(^_Jy3+v+mF4x0ZQscOD|c;eQ+rVNh|!*$?G zJSJe{^J`yjB3?qEi}rx$aJgW2v4T%&oP6%awz}8+aI8RBQ~Gs!b{B!K*0ZG)nIuSn z$`%Kc?z!05$oh3YLMNb><_!5-*T>+H5UH=&Y1GUsrxs#)SjQbM4w%^fB+Lg_D1e__ z0rc9Mo+%~5jp{w0*i(l^alun#ezO1Q^;Weh*+ZC{Hxe3eE+{&K$jg~(2kZ^Q3;)JJ z{|E|79?j;E7o${e-CHpR6+xaS zo)SKZGZmUPl@9yU>$3CF_(nAYFU^U{O8RvpwNUwgqv+Wun#C5!ix~m+lU;*CIQ=IA z3z_3V#Sg$+xD&vKTpl#fu0hY?we@3TR8>uc?(wXP7a3AixPR2lY}*NdQrJ-xAb)~O z6Uk;D;kK_)eFn5ZY{h-wuc&I*$25>``TDv?>`rr8Ci;B*mTQ`2HQD~pR{58OFMUggU$DIyaVxQ1>lZQW zrDLUy4;|6#9f}|!9>wEn<`%f-4UtxjKOHVmV1R%{Y)HR68FoI+4Mcv7kwsDI-jPYR z2ahO(PWNxTLGA_j7)`z;KNZD=`6fK) zO@*LL`fS~AOY=u_lTmfoY0&)SqH;f&{kU#Y^V`g*Q0FP$Gh{t*r%H1c>7{ ze~I&6e|^=HQ$u8M4SRo$<{3~R+n`8Sl~vB#ah~FXW@!|vl~RN_YefD%<{VJO7gMR` zK~RawnStGCk67!qf!q@xn^Dm!$gkmj`!@?jDR0$Hc6V!4sE#dkH$OUDF>)oPA&&p4 zg}FZnG;(Uy)ZREfR*r06cCfv!c2;$T=UF=CzStDL>@{@6h`a5PA3c9s`@VOJl`ffA z#2#WapiODgl401EhCd7jj9X9iM!2NdN?G21-Ybj8M1pTN`LxJScc->lY=5o&p~~>Km|0+w&nEi?5wEBWXZz$r)9#0Er=RR00w}d z*34sVn5%uP(T3R6`n}x({PrL6} zqwqSySLnkRrc9dT74+7&S_fyUOrX1>A>(w}M%xSZW<#B>i-T9!6fGU!JnN;#l|oB{ z9lJFKleHbUnVL?Rh{o7eUhmAZzwABcmQ z^tTTnKN3`f1V5#(qVn#*ffm+J;DWUum}hSO#&o0WpUuV&q2GsWbiH7#SbcEtDKe)u z_QMVAS<2N7MklF7sI9r|^+uw(=@ZCYMXw~u+LHMX1u=B9xa;v^7`9AgJ(|@o!!C~4 z==(IvUvEu)#OuE4znH%%XmpGjDUPKvK`N$h$|alB*-DLHzgBO0M)-5fa-5pltBw^y zFqcg)4->6sbT1 z{sOvto3{`jz2fYVs{kLoY^aAO@R;%kMOi8PAcZ7&&_g3Z&=%!1PfSXcZ#hr6fQn>(%@p$hL+ zwqL*GmaTTL?#>MYINt51RXyXeO^z$QeNqHU4@ej2$UVREzR`vNUzn934vkvXv#OQV zGYYC_Jxq$^QLw@(zJ zbzd6#obtl^3I4XexFGCJ_FM3Sc07vdz7fbpt{9hi`0#V$lA79&^9XRxe>%_}VW%-5 zTo8u<)X7Rhko+KBWE$_EhL)7fMqGNW{ib;&$$-$h^KYoGGe5N@Dj;a|=~f})rB{W7 za7>kh4Q3zr47UDR!?kEzp}CDs(C8YI$AiYX)KNYwynX~1hL;>XcyM~7@Q8sum=I>) z;nkPxUc2S_LV94#wFL!xOP~yggvTCB*LE&k-V&)%X~Dl8rD{PLie?N+aPkdkGzvQ+6(tWCJal zntmDIPAtJE$c7&y(8e{`w*2kyFI=AHKvf*Rs;aaKdiO-(t-Y4^Fr-0+&^UX1 z%IO>dx8$9?uJ4Fze}f(OPseKr|DhppGXH-AB>&gQ%KxJ`9v}KqB(m~){`@&%j{A`c zx{*OxBi?oSgf_i2`sxZYy}F}ZTk(`#%!CD!O_UhZ8dn!B^`cMESn=%n^Nf~RZD<(W zHlJV_zSzoflY&i$K#_trl)M%F;mvRwZ>OUfy-%y`c97Ox`@3@6m)t=l_Tdldw$G~BoQoy&-P zfX##na5tF8vf(`qM+=OI&%HtmAp&E~@`E}nY5l;cu?qHo9yT!+xK4zD#BDvg&Y)65 zum%Q&l&3E!3>`znOkwdqaAEY-&^uzT6om@Zo)OcWTcVQ5PU)x$fJ#VMg-nAmof-TR zYXoUdU|KLF$=<+Y#*AyfQ7I{o9FC-!z~HRoRBnYZcKvz|poN1*7?^T%Cd~3m~nwNY4SH1M?MiK}7j_t2s-RdOOMPFub-1uj+ zJmODE+$p{q@NSSi);twCDd~-27dm~wA*>rmpr!v0G z2gMqlLg)Ce)>aMI;jgQ!N6AT!QJBa!7&-DPd8eK*h_Yfk;Mn|=$6r6q!A2ou;w@P4 z)i~3}+PAoc>`4_%_93K_(U!}u8Yv?T8_oHFKc?@Ztm=j?xQ(QpxSp@6XakR_IbBiX zl+^|ziLEx)g%@?!(dEHr8VBM8x|3?Nbn-c28|>qb0-$Re>(c2&iT=g&=Xr?6kX_QA z!oIG;{S4Jxll~>DcYhU>jD5)AAa62HdJEf)%!$HdKA8)+GZ-|pFpjdRE}lELd*42^ z92Nj{()$Bq&EOEv&`a^kJi}3@r_S@ z9FiuSV{iDk`5mFy0!bO ziBF{9D%h$Qd!mTqN@Y(N0Fl_XkumV8siS3(<8%%@GDXeEnB5emAVj|PJxfTO-|A!c z*PpPH+i-38^^wr?Ox1?zePOV1GvynwKre$kG{6PFU5$`ojX4_#><$ZcV7+EnD?LJffxvT?F%HckJNP<|F87cM_RW0< z6+HD}dX5*@3~7EQ=%hI4czAhz1uX!!%)C+sAp7FkGqpb64y9qnDL}&~=bkj>G?85o z)o}|e+*qVB_6z~>lIuIXs^sA6!cQ~H9~ImnP*t?D|JlmPNo)IxoQs2^hNh-4Z@cg~ z30`TE-fIL)Wvq$5qR$Mi*6ST=R3nNYHDiEcA${UQ{KjOpb#u@dhSq|Ke}Qs!rtzp0 z^ii@3#7hC4vzJYm)R@1@za0T*JV5n`;PYhj@U|QprMmbiWI{3;RB^zYOgX>x!`VQ04g+o5s0w=Z1w%2}$uj;HWvj5+3uim}}X`XZ`NQTiYXimr2aRzfx^z z^rAN4=H`YUfQAN&iN05kh0UGmH}lv=8cAmy&zjvj$AVPxd*`gv$f#1W_KK4YVBUUOpb+y zn)Wr;HL&!~ofxU1^789XAKDNdeR^rk!KZNHyqMI4OB4$BB>aF9$bCrm>2pK*=GA4C zS-$nSq$%hg+gCoiV{H(0EWle_Xx4fJAUs+*dW&^Py7FP|C51d;f=w)}0S5~;XsXS{ zB_(R3M$M?XpCO5p9{1V&(SlX#pRZof0l<02ZiT@JSE>fxD^90b{)FJW5Lc45C`hs_ z6gGs=SIKZrtrPl1;j>5fD|zB)x9kyExtCuzftY|Fj;g{G*p=!#jLeP6s6XFn;Y6P; z;Tf4fBIl3cOJW*4RuqVwcWgT#YE6MIYAi=o7<$XX4mi@@Uzs@&+A~+;7+4lL{ z7noE-Xf7+(%eGIaD-;-=ANMy9282G4A)`x|eYfyhzy8g*{7D>-ppGU5i~6p?V=2b9 z7NX99^AQRd8X6${)F2YNystPb5Y_p&B=k5voVm^@mf^-+3b3T0v)G1t%IcdqDmf(l)cv!wjiz*d~nLZf2mY z)Tds-y}QkK7}`~|pu$9n;xm~Ffv5WWYF&@<+z6rI-eHp(ImyP+2?_0B)7V+?4sY{m zD@S?I0mXPHK?y##=AzK*NGu16PgO0Ka*E`>T~tT+(!fpoMOGG}wX)w!uNwR02j9N= zut60`Z^e$=((1Wcmsx7E;g1@3h+gdcb-l~CXyokQ&3ds7xU9Zx#KZ?3`aOcr2u0m52zhQx$);(9Y^Nt7=<%YgJNviU%y|?Grf2ZpFuCvRm XggXXx6X&hqouWA=^Ncf%7H|C@+bX0o From 935f40c81c97fbeb38a43266a1c87aeb22979316 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:51:33 -0300 Subject: [PATCH 05/29] feat(linalg): Trying using MINPACK --- src/linalg.f90 | 46 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/src/linalg.f90 b/src/linalg.f90 index cda3bf2..0e9d2f8 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -185,6 +185,7 @@ subroutine full_newton(fun, iters, X, ns, S, max_iters, F, dF) !! !! Procedure that solves a point with the Newton-Raphson method. use constants, only: ouput_path + use minpack_module, only: dpmpar, hybrj1 interface subroutine fun(X, ns, S, F, dF) !! Function to solve @@ -207,28 +208,59 @@ subroutine fun(X, ns, S, F, dF) !&> real(pr) :: b(size(X)), A(size(X), size(X)) - real(pr) :: dX(size(X)), tol = 1e-2 - integer :: funit_log_newton + real(pr) :: dX(size(X)), tol = 1e-6 - integer :: n + integer :: n, ldfjac, info, lwa, funit + real(pr) :: fvec(size(x)), fjac(size(x),size(x)) + real(pr), allocatable :: wa(:) n = size(X) + ! ldfjac = n + ! lwa = (n*(n+13))/2 + ! allocate(wa(lwa)) + ! call hybrj1(fcn, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa) + + ! f = fvec + ! df = fjac + ! iters = 3 + ! if (info == 4) iters = max_iters+1 + dX = 20 + ! open(newunit=funit, file="fenvelopes_output/newton") + ! write(funit, *) n + b = 500 + newton: do iters = 1, max_iters - if (maxval(abs(dx/x)) < tol) exit newton + if (maxval(abs(dx/x)) < tol .or. maxval(abs(b)) < tol) exit newton call fun(X, ns, S, b, a) b = -b + + if (any(isnan(b))) dx = dx/maxval(abs(dx)) dX = solve_system(A, b) - do while (maxval(abs(dx)) > 0.5*maxval(abs(x))) - dX = dX/2 + do while(maxval(abs(dx)) > 0.5) + dX = dX/10 end do + ! write(funit, *) x, dx, -b + X = X + dX end do newton F = -b dF = A + ! close(funit) + + contains + subroutine fcn(n, x, fvec, fjac, ldfjac, iflag) + integer, intent(in) :: N + real(pr), intent(in) :: x(n) + real(pr), intent(inout) :: fvec(n) + real(pr), intent(inout) :: fjac(ldfjac, n) + integer, intent(in) :: ldfjac + integer, intent(inout) :: iflag + call fun(X, ns, S, fvec, fjac) + end subroutine end subroutine -end module linalg +end module linalg \ No newline at end of file From 12225dc4670bb9b3072a7de1b44f02aa2f15b84b Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:51:45 -0300 Subject: [PATCH 06/29] form --- src/linalg.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg.f90 b/src/linalg.f90 index 0e9d2f8..ff1a49c 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -104,7 +104,7 @@ function intersect_one_line(lx, ly) result(intersections) allocate (intersections(0)) line1: do i = 2, size(lx) - 1 - line2: do j = i + 15, size(lx) + line2: do j = i + 2, size(lx) associate ( & x1 => lx(i - 1), x2 => lx(i), & x3 => lx(j), x4 => lx(j - 1), & From 066e9a3411af02dd8d5271f6c9522d612eb0510f Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:52:08 -0300 Subject: [PATCH 07/29] feat(fpm): added minpack dependency --- fpm.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/fpm.toml b/fpm.toml index a71612a..9c6786f 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 = ["lapack", "minpack"] +link = ["lapack"] auto-executables = true auto-tests = true auto-examples = true @@ -16,6 +16,7 @@ library = false [dependencies] stdlib = "*" +minpack= "*" ftools = { git = "https://github.com/fedebenelli/ftools" } yaeos = { git = "https://github.com/fedebenelli/yaeos" } FLAP = { git = "https://github.com/szaghi/FLAP", tag="v1.2.5" } From f88e7717dfb809f28233e046de79a03b2ceaeca9 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:53:13 -0300 Subject: [PATCH 08/29] feat(Kinit): Kwilson for bubble points Added the possibility to specify t0 and p_end on Kwilson calculation --- src/new/envelopes.f90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index 384b0a4..6b771ec 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -35,20 +35,21 @@ function F(X, ns, S) ! =========================================================================== ! Initializators ! --------------------------------------------------------------------------- - subroutine k_wilson_bubble(z, t, p, k) + subroutine k_wilson_bubble(z, t_0, p_end, t, p, k) !! Find the Wilson Kfactors at ~10 bar to initialize a bubble point - ! use system, only: pc, tc, w use legacy_ar_models, only: pc, tc, w real(pr), intent(in) :: z(:) + real(pr), intent(in) :: t_0 + real(pr), intent(in) :: p_end real(pr), intent(in out) :: p real(pr), intent(in out) :: t real(pr), intent(out) :: k(size(z)) P = 100.0 - T = 200.0 + T = t_0 - do while (P > 10) + do while (P > p_end) T = T - 5._pr P = 1.0_pr/sum(z * pc*exp(5.373_pr*(1 + w)*(1 - tc/T))) end do From 89c3d3c564c047920e7b89fa411ecb1f264fdbe9 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:54:01 -0300 Subject: [PATCH 09/29] feat(F2PT): Added default cases Calculate stable phase when no incipient phase is given --- src/new/envelopes.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index 6b771ec..47e22da 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -163,7 +163,7 @@ function X2(kfact, P, T) result(X) end function subroutine F2(incipient, z, y, X, S, ns, F, dF) - character(len=:), allocatable, intent(in) :: incipient + character(len=*), intent(in) :: incipient real(pr), intent(in) :: z(:) real(pr), intent(in) :: X(nc + 2) real(pr), intent(in) :: y(nc) @@ -199,6 +199,9 @@ subroutine F2(incipient, z, y, X, S, ns, F, dF) case ("2ndliquid") ix = 1 iy = 1 + case default + ix = 0 + iy = 0 end select call TERMO(n, iy, 4, T, P, y, Vy, lnfug_y, dlnphi_dp_y, dlnphi_dt_y, dlnphi_dn_y) From 511f7184194e9303c83c4053e0cd6715a9a0d802 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:56:38 -0300 Subject: [PATCH 10/29] feat(HPLL): Generalizing method to find HPLL line --- src/new/envelopes.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index 47e22da..dfcff05 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -262,21 +262,23 @@ subroutine find_hpl(t, p, k) real(pr), intent(in) :: p real(pr), intent(out) :: k(nc) - integer :: i + integer :: i, ncomp real(pr) :: diff real(pr) :: v real(pr) :: x(nc), y(nc), lnfug_z(nc), lnfug_y(nc) diff = -1 + ncomp = nc + y = 0 - y(nc) = 1 + y(ncomp) = 1 do while(diff < 0) t = t - 1.0_pr call termo(nc, 4, 1, t, p, z, v, philog=lnfug_z) call termo(nc, 4, 1, t, p, y, v, philog=lnfug_y) - diff = (log(z(nc)) + lnfug_z(nc)) - (log(y(nc)) + lnfug_y(nc)) + diff = (log(z(ncomp)) + lnfug_z(ncomp)) - (log(y(ncomp)) + lnfug_y(ncomp)) end do k = exp(lnfug_y - lnfug_z) From a9612fee04210ef6c2469f5269d09f32a4acfa12 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 09:58:01 -0300 Subject: [PATCH 11/29] refactor(PT2ph): Trying different limits --- src/new/envelopes.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index dfcff05..bb37691 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -619,21 +619,21 @@ subroutine envelope2(ichoice, n, z, T, P, KFACT, & ! This will probably always e T = exp(X(n + 1)) - if (.not. passingcri .and. abs(T - Told) > 7) then + do while (.not. passingcri .and. abs(T - Told) > 15) ! Delta T estimations > 7K are not allowed delS = delS/2 S = S - delS X = Xold + dXdS*delS T = exp(X(n + 1)) - end if + end do P = exp(X(n + 2)) KFACT = exp(X(:n)) y = z*KFACT ! Finish conditions - if ((dXdS(n + 1)*delS < 0 .and. P < 0.1 .and. T < 120.0) & ! dew line stops when P<0.1 bar or T<150K - .or. (P > 1.0 .and. T < 50.0) & ! bubble line stops when T<150K + if ((dXdS(n + 1)*delS < 0 .and. P < 0.05 .and. T < 50.0) & ! dew line stops when P<0.1 bar or T<50K + .or. (P > 1.0 .and. T < 50.0) & ! bubble line stops when T<50K .or. (P > 5000) & .or. (abs(dels) < 1.d-10)) then run = .false. From 2b30995c950ab8fd7c662f8bafc896fe238081ca Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:01:28 -0300 Subject: [PATCH 12/29] refactor(PT2ph): Detections - Avoiding stepping back on CP detection due to sometimes getting stuck there - Detection of new crossing cases --- src/new/envelopes.f90 | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index bb37691..fe9d86e 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -604,16 +604,18 @@ subroutine envelope2(ichoice, n, z, T, P, KFACT, & ! This will probably always e ! the step given by the most changing logK to fall into ! the black hole - passingcri = .true. - if (stepX > 0.07) then - ! half step back - S = S - delS/2 - X = X - dXdS*delS/2 - else - ! one more step to jump over the critical point - S = S + delS - X = X + dXdS*delS - end if + S = S + delS + X = X + dXdS*delS + ! passingcri = .true. + ! if (stepX > 0.07) then + ! ! half step back + ! S = S - delS/2 + ! X = X - dXdS*delS/2 + ! else + ! ! one more step to jump over the critical point + ! S = S + delS + ! X = X + dXdS*delS + ! end if end do end block critical_region @@ -1035,12 +1037,15 @@ subroutine get_case(& else if (size(inter_hpl_bub) == 1 .and. size(inter_dew_bub) == 1) then this_case = "2_HPL_BUB_DEW_BUB" intersections = [inter_hpl_bub, inter_dew_bub] - else if (size(inter_hpl_bub) == 1) then - this_case = "1_HPL_BUB" - intersections = inter_hpl_bub + else if (size(inter_hpl_bub) == 2) then + this_case = "2_HPL_BUB" + intersections = [inter_hpl_bub(1), inter_hpl_bub(2)] else if (size(inter_hpl_dew) == 1) then this_case = "1_HPL_DEW" intersections = inter_hpl_dew + else if (size(inter_hpl_bub) == 1) then + this_case = "1_HPL_BUB" + intersections = inter_hpl_bub else this_case = "0" allocate(intersections(0)) From 6c271a28406f93c3678ab313ef2dbce51c828da2 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:02:08 -0300 Subject: [PATCH 13/29] feat(PT3ph): Break on positive or negative beta --- src/new/envelopes.f90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index fe9d86e..606f80f 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -996,15 +996,17 @@ function break_conditions_three_phases(X, ns, S) real(pr) :: S !! Value of specification integer :: n - real(pr) :: p, t + real(pr) :: p, t, beta logical, allocatable :: break_conditions_three_phases(:) n = (size(X) - 3)/2 p = exp(X(2*n + 1)) + beta = x(2*n+3) t = X(2*n + 2) - - break_conditions_three_phases = [ .false.& - ! p < 1 .or. p > 5000 & + + break_conditions_three_phases = [ & + beta < 0 & + .or. 1 < beta & ] end function ! =========================================================================== @@ -1106,7 +1108,7 @@ subroutine pt_three_phase_from_intersection(& X = [lnKx, lnKy, log(p), log(t), beta] call pt_envelope_three_phase(X, ns, del_S0, pt_x_3(i_inter)) ! ================================================================== - + ! ================================================================== ! Line with incipient phase liquid ! ------------------------------------------------------------------ From f12e5f40ef6b1696327116e34cf631b6be7fe2a7 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:02:41 -0300 Subject: [PATCH 14/29] feat(PT3ph): Critical points factors Adjusted CP detection and jump factors, should be revisited --- src/new/envelopes.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index 606f80f..9282580 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -927,7 +927,7 @@ subroutine pt_envelope_three_phase(X0, spec_number, del_S0, envel) Xnew(size(X0)), fact real(pr) :: pc, tc, dS_c, dXdS_in(size(X0)) integer :: max_changing, i - fact = 3.0_pr + fact = 4.0_pr loop: do i = 0, 1 Xnew = X + fact*dXdS*del_S @@ -948,8 +948,7 @@ subroutine pt_envelope_three_phase(X0, spec_number, del_S0, envel) pc = exp(Xnew(2*n + 1)) cps = [cps, critical_point(tc, pc, 0.0_pr)] - del_S = dS_c + sign(0.04_pr, dS_c) ! dS_c + 2*del_S ! * fact - ! del_S = del_S * fact + del_S = dS_c + sign(0.04_pr, dS_c) write (funit_output, *) "" write (funit_output, *) "" @@ -958,7 +957,6 @@ subroutine pt_envelope_three_phase(X0, spec_number, del_S0, envel) end do loop end block detect_critical - if (x(2*n + 3) > 1 .or. (x(2*n+3) < 0)) exit enveloop X = X + dXdS*del_S S = X(ns) From c9db313ef7c010d19546c70a21cae241ecef3d23 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:03:04 -0300 Subject: [PATCH 15/29] feat(PT2ph): Increasing deltaS jump --- src/new/envelopes.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/new/envelopes.f90 b/src/new/envelopes.f90 index 9282580..2086a01 100644 --- a/src/new/envelopes.f90 +++ b/src/new/envelopes.f90 @@ -124,7 +124,7 @@ subroutine update_specification(iter, passingcri, X, dF, ns, S, delS, dXdS) else delS = max(updel, -delmax) end if - delS = delS*3 + delS = 5*delS S = S + delS end subroutine From 8fdbe61872c1624cb005e6fd89040fcdd6d416d7 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:09:27 -0300 Subject: [PATCH 16/29] prints on output after namelist reading --- src/new/io_nml.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/new/io_nml.f90 b/src/new/io_nml.f90 index 9f8d762..2e161c4 100644 --- a/src/new/io_nml.f90 +++ b/src/new/io_nml.f90 @@ -260,9 +260,9 @@ subroutine write_system(file_unit) write (file_unit, *) "====================" write (file_unit, *) "EoS Parameters: " write (file_unit, *) "--------------------" - write (file_unit, fmt_pure) "ac: ", ac - write (file_unit, fmt_pure) "b: ", b - write (file_unit, fmt_pure) "k: ", k + write (file_unit, fmt_pure) "ac:", ac + write (file_unit, fmt_pure) "b :", b + write (file_unit, fmt_pure) "k :", k write (file_unit, *) "====================" write(file_unit, *), "Mixing Rules: " @@ -272,12 +272,12 @@ subroutine write_system(file_unit) write(file_unit, fmt_names) names do i = 1, nc - write (file_unit, fmt_pure) trim(names(i)), kij(i, :) + write (file_unit, fmt_pure) str(i), kij(i, :) end do write(file_unit, *) "lij: " write(file_unit, fmt_names) names do i = 1, nc - write (file_unit, fmt_pure) trim(names(i)), lij(i, :) + write (file_unit, fmt_pure) str(i), lij(i, :) end do end subroutine end module From 548ce6e829931645109ec1b6cbafed08ce0c9170 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:12:39 -0300 Subject: [PATCH 17/29] WIP(Px2ph): Initialization from PT --- src/new/mod_inj_envelopes.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index cab7c77..e49eb11 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -779,7 +779,7 @@ function px_two_phase_from_pt(t_inj, pt_env_2, t_tol, del_S0) result(envel) real(pr) :: del_S del_S = optval(del_S0, 0.1_pr) - pold = 0 + pold = 1e9 ts_envel = pack(pt_env_2%t, mask=abs(pt_env_2%t - t_inj) < t_tol) do i = 1, size(ts_envel) @@ -789,7 +789,7 @@ function px_two_phase_from_pt(t_inj, pt_env_2, t_tol, del_S0) result(envel) pt_env_2%p(idx), pt_env_2%p(idx + 1), & t_inj) - if (abs(p - pold) < 5) cycle + if (abs(p - pold) < 2) cycle pold = p k = exp(interpol( & From f650b02c97130bc0474dbb74be156bc6cfa13b43 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:16:54 -0300 Subject: [PATCH 18/29] feat(Px3ph): Added more details to ouput Also pointer variables for readability --- src/new/mod_inj_envelopes.f90 | 29 +++++++++++++++++++++++------ 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index e49eb11..bc0ee69 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -308,9 +308,9 @@ subroutine fix_step_two_phases(X, ns, S, solve_its, del_S, dXdS) real(pr) :: dP, dalpha real(pr) :: dP_tol, dalpha_tol integer :: n - + n = size(X) - 2 - + del_S = sign(del_S_multiplier, del_S)*minval([ & max(sqrt(abs(X(ns)))/10, 0.1), & abs(del_S)*3/solve_its & @@ -486,10 +486,12 @@ subroutine F_injection_three_phases(Xvars, ns, S, F, dF) ) df(i + n, 2*n + 2) = sum(Ky*dlnphi_dn_y(i, :)*dwda & - dlnphi_dn_w(i, :)*dwda) + ! Derivatives wrt beta df(i, 2*n + 3) = sum(Kx*dlnphi_dn_x(i, :)*dwdb & - dlnphi_dn_w(i, :)*dwdb) df(i + n, 2*n + 3) = sum(Ky*dlnphi_dn_y(i, :)*dwdb & - dlnphi_dn_w(i, :)*dwdb) + df(2*n + 1, i) = Kx(i)*dwdKx(i) df(2*n + 1, i + n) = Ky(i)*dwdKy(i) @@ -526,7 +528,7 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) type(critical_point), allocatable :: cps(:) - real(pr) :: X(size(X0)) + real(pr), target :: X(size(X0)) integer :: ns real(pr) :: S real(pr) :: XS(max_points, size(X0)) @@ -534,6 +536,13 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) real(pr) :: F(size(X0)), dF(size(X0), size(X0)), dXdS(size(X0)) + real(pr), pointer :: lnP + real(pr), pointer :: alpha + real(pr), pointer :: beta + real(pr), pointer :: lnKx(:) + real(pr), pointer :: lnKy(:) + real(pr) :: z(size(z_0)) + integer :: point, iters, n integer :: i integer :: funit_output @@ -546,6 +555,13 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) S = X(ns) del_S = del_S0 + lnKx => X(:n) + lnKy => X(n+1:2*n) + lnP => X(2*n+1) + alpha => X(2*n+2) + beta => X(2*n+3) + call get_z(alpha, z) + ! ====================================================================== ! Output file ! ---------------------------------------------------------------------- @@ -558,7 +574,7 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) open (newunit=funit_output, file=fname_env) write (funit_output, *) "#", T write (funit_output, *) "STAT", " iters", " ns", " alpha", " P", & - " beta", (" lnKx"//str(i), i=1,n), (" lnKy"//str(i), i=1,n) + " beta", (" lnKx"//str(i), i=1,n), (" lnKy"//str(i), i=1,n), (" z" //str(i), i=1,n) write (funit_output, *) "X0", iters, ns, X(2*n + 2), exp(X(2*n + 1)), & X(2*n + 3), X(:2*n) ! ====================================================================== @@ -572,8 +588,9 @@ subroutine injection_envelope_three_phase(X0, spec_number, del_S0, envels) exit enveloop end if - write (funit_output, *) "SOL", iters, ns, X(2*n + 2), & - exp(X(2*n + 1)), X(2*n + 3), X(:2*n) + call get_z(alpha, z) + write (funit_output, *) "SOL", iters, ns, alpha, & + 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) From d932deb24a16d1833af573085703499f0a937b09 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:18:38 -0300 Subject: [PATCH 19/29] stepping parameters --- src/new/mod_inj_envelopes.f90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index bc0ee69..11b16a3 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -17,22 +17,23 @@ module inj_envelopes ! Parameters ! --------------------------------------------------------------------------- integer :: env_number = 0 !! Number of calculated envelope - integer :: max_iters = 150 !! Maximum number of iterations for a newton step + integer :: max_iters = 50 !! Maximum number of iterations for a newton step integer, parameter :: max_points = 1000 !! Maximum number of points for each envelope real(pr), allocatable :: z_0(:) !! Original fluid composition real(pr), allocatable :: z_injection(:) !! Injection fluid composition real(pr) :: T !! Temperature of injection character(len=10) :: injection_case !! Kind of injection displace|dilute character(len=255) :: FE_LOG - + ! Two-phase settings - real(pr) :: del_S_multiplier = 1.6_pr + real(pr) :: del_S_multiplier = 2.0_pr real(pr) :: max_dp=50.0_pr real(pr) :: max_dalpha=0.01_pr + real(pr) :: critical_multiplier = 2.0_pr ! Three phase parameters - real(pr) :: del_S_multiplier_three_phase = 1.5_pr - real(pr) :: critical_fact = 4.0_pr + real(pr) :: del_S_multiplier_three_phase = 1.7_pr + real(pr) :: critical_fact = 3.0_pr ! =========================================================================== contains subroutine from_nml(filepath) From 18d407825abf06230af55a7cfa0ee7691107dfca Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:20:25 -0300 Subject: [PATCH 20/29] critical points --- src/new/mod_inj_envelopes.f90 | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index 11b16a3..1c4b880 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -209,7 +209,7 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) Xnew(size(X0)), fact real(pr) :: pc, alpha_c, dS_c integer :: max_changing - fact = 2.5 + fact = critical_fact ! 2.5 Xnew = X + fact*dXdS*del_S @@ -220,16 +220,17 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) max_changing = maxloc(abs(K - Knew), dim=1) dS_c = ( & - -k(max_changing)*(Xnew(ns) - X(ns)) & + -K(max_changing)*(Xnew(max_changing) - X(max_changing)) & /(Knew(max_changing) - K(max_changing)) & ) - del_S = dS_c*1.1 - Xnew = X + dXdS*dS_c + Xnew = X + dXdS * dS_c alpha_c = Xnew(n + 2) - pc = Xnew(n + 1) - + 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 @@ -253,7 +254,7 @@ subroutine injection_envelope(X0, spec_number, del_S0, envels) write (funit_output, *) "#critical" if (size(cps) > 0) then do i = 1, size(cps) - write (funit_output, *) cps(i)%t, cps(i)%p + write (funit_output, *) cps(i)%alpha, cps(i)%p end do else write (funit_output, *) "NaN NaN" From 163302536f2df5ddba5df9fe72a6c2619536f991 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:20:58 -0300 Subject: [PATCH 21/29] initializatoin from DSP --- src/new/mod_inj_envelopes.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index 1c4b880..e6cd175 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -280,6 +280,7 @@ subroutine update_spec(X, ns, del_S, dF, dXdS) real(pr) :: dFdS(size(X)) integer :: ns_new + integer :: n dFdS = 0 dFdS(size(dFdS)) = 1 @@ -294,7 +295,6 @@ subroutine update_spec(X, ns, del_S, dF, dXdS) dXdS = dXdS/dXdS(ns_new) ns = ns_new end if - end subroutine subroutine fix_step_two_phases(X, ns, S, solve_its, del_S, dXdS) @@ -915,7 +915,7 @@ function px_three_phase_from_inter(& integer :: ns del_S = optval(del_S0, -0.05_pr) - beta = optval(beta0, 1.0_pr) + beta = optval(beta0, 0.99_pr) i = inter%i j = inter%j From 5bff29eb99eb70b5a4dc40835c918b27c2b116cd Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:22:20 -0300 Subject: [PATCH 22/29] limiting a bit the calculation of isolated lines, tihs should be perfected later --- src/new/mod_inj_envelopes.f90 | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/src/new/mod_inj_envelopes.f90 b/src/new/mod_inj_envelopes.f90 index e6cd175..9954e9b 100644 --- a/src/new/mod_inj_envelopes.f90 +++ b/src/new/mod_inj_envelopes.f90 @@ -986,14 +986,17 @@ function px_hpl_line(alpha_0, p) find_hpl: do ncomp = nc, 1, -1 alpha_in = alpha_0 p_in = p - + y = 0 y(ncomp) = 1.0 - - diff = foo([alpha_in, p_in, y] ) - do while (abs(diff) > 0.01 .and. p_in < 5000) + + diff = foo([alpha_in, p_in, y]) + return + ! print *, ncomp, diff + do while (abs(diff) > 0.001 .and. p_in < 5000) p_in = p_in + 10.0_pr diff = foo([alpha_in, p_in, y]) + ! print *, diff end do if (p_in >= 5000) then @@ -1003,8 +1006,12 @@ function px_hpl_line(alpha_0, p) end do end if - if (alpha_in > 0 .or. p_in < 5000) then - call optim + p_in = 600 + alpha_in = 0.9 + + 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) X = [log(K), log(P_in), alpha_in] @@ -1022,7 +1029,7 @@ function foo(x) result(f) real(pr) :: f if (x(1) > 1) x(1) = 0.97_pr - alpha_in = x(1) + ! alpha_in = x(1) p_in = x(2) y = abs(x(3:)) From d8aced91b2dece78809494d55034469c5e003ce4 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:23:53 -0300 Subject: [PATCH 23/29] Don't remove lowP PT; this can give useful lines too! --- app/main.f90 | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index 266f17f..3efade4 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -131,19 +131,6 @@ subroutine pt_envelopes n_points, Tv, Pv, Dv, ncri, icri, Tcri, Pcri, Dcri, & pt_dew & ) - - ! Remove the low pressure parts. - ! n = 1 - ! do i = 2, size(pt_dew%t) - ! n = n + 1 - ! if (pt_dew%t(i) - pt_dew%t(i - 1) < 0) exit - ! end do - - ! if (n /= size(pt_dew%t)) then - ! pt_dew%t = pt_dew%t(i:) - ! pt_dew%p = pt_dew%p(i:) - ! pt_dew%logk = pt_dew%logk(i:, :) - ! end if ! ======================================================================== ! ======================================================================== From c66fcaaa6201d0a6b80d671dd602184ff02aca9d Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:24:17 -0300 Subject: [PATCH 24/29] Kwilson initialization --- app/main.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index 3efade4..617d8cc 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -106,12 +106,12 @@ subroutine pt_envelopes ! ======================================================================== ! Bubble envel ! ------------------------------------------------------------------------ - call k_wilson_bubble(z, t, p, k) + call k_wilson_bubble(z, t_0=230.0_pr, p_end=0.5_pr, t=t, p=p, k=k) call envelope2( & 1, nc, z, T, P, k, & n_points, Tv, Pv, Dv, ncri, icri, Tcri, Pcri, Dcri, & pt_bub & - ) + ) ! ======================================================================== ! ======================================================================== From 5c189b1ba4590b8d5ed026784e7703b4fb889065 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:24:42 -0300 Subject: [PATCH 25/29] Dew line initial temperature. This should be remove later --- app/main.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/main.f90 b/app/main.f90 index 617d8cc..aaa884e 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -117,7 +117,7 @@ subroutine pt_envelopes ! ======================================================================== ! Dew/AOP envelopes ! ------------------------------------------------------------------------ - t = 200 + t = 300 p = p_wilson(z, t) do while (p > 0.1) t = t - 5 From 04e4365e4beb089bd2d129f08a3bd238f38eda6a Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:25:34 -0300 Subject: [PATCH 26/29] Possibility of calculate PT line on a specific alpha. Should be refacetored! --- app/main.f90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index aaa884e..34a67b8 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -1,9 +1,9 @@ program main - use envelopes, only: PTEnvel3 use dtypes, only: envelope - use inj_envelopes, only: injelope + use envelopes, only: PTEnvel3 + use inj_envelopes, only: injelope, get_z use constants, only: pr, ouput_path - use legacy_ar_models, only: nc + use legacy_ar_models, only: nc, z use flap, only: command_line_interface use stdlib_ansi, only: blue => fg_color_blue, red => fg_color_red, & operator(//), operator(+), & @@ -20,9 +20,12 @@ program main type(PTEnvel3), allocatable :: pt_bub_3(:), pt_dew_3(:) !! Shared 3ph-PT envelopes type(injelope) :: px_bub, px_dew, px_hpl !! Shared 2ph-Px envelopes + real(pr) :: alpha=0.0 + ! Setup everything call setup + call get_z(alpha, z) ! PT Envelopes call cpu_time(st) call pt_envelopes From dff2b5eb0f20016189c5dadae8f471ef3befe8d2 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:25:48 -0300 Subject: [PATCH 27/29] Calculation of DSP lines --- app/main.f90 | 42 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/app/main.f90 b/app/main.f90 index 34a67b8..eef673b 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -166,10 +166,17 @@ subroutine pt_envelopes three_phase: block use envelopes, only: pt_three_phase_from_intersection - allocate(pt_bub_3(size(intersections)), pt_dew_3(size(intersections))) select case(pt_case) case("2_DEW_BUB") + dsp_line: block + use dsp_lines, only: injelope, dsp_line_from_dsp + type(injelope):: dsps(2) + integer :: i + do i=1,size(intersections) + dsps = dsp_line_from_dsp(intersections(i), pt_dew, pt_bub) + end do + end block dsp_line call pt_three_phase_from_intersection(& pt_dew, pt_bub, intersections, & pt_bub_3, pt_dew_3 & @@ -183,12 +190,45 @@ subroutine pt_envelopes pt_dew, pt_bub, [intersections(2)], & pt_bub_3, pt_dew_3 & ) + case("2_HPL_BUB") + call pt_three_phase_from_intersection(& + pt_hpl, pt_bub, [intersections(1)], & + pt_bub_3, pt_dew_3 & + ) + call pt_three_phase_from_intersection(& + pt_hpl, pt_bub, [intersections(2)], & + pt_bub_3, pt_dew_3 & + ) + dsp_line_2hpl_bub: block + use dsp_lines, only: injelope, dsp_line_from_dsp + type(injelope):: dsps(2) + integer :: i + do i=1,size(intersections) + dsps = dsp_line_from_dsp(intersections(i), pt_hpl, pt_bub, alpha0=alpha) + end do + end block dsp_line_2hpl_bub case("1_HPL_DEW") + dsp_line_hpl: block + use dsp_lines, only: injelope, dsp_line_from_dsp + type(injelope):: dsps(2) + integer :: i + do i=1,size(intersections) + dsps = dsp_line_from_dsp(intersections(i), pt_hpl, pt_dew, alpha0=alpha) + end do + end block dsp_line_hpl call pt_three_phase_from_intersection(& pt_hpl, pt_dew, intersections, & pt_bub_3, pt_dew_3 & ) case("1_HPL_BUB") + dsp_line_hpl_bub: block + use dsp_lines, only: injelope, dsp_line_from_dsp + type(injelope):: dsps(2) + integer :: i + do i=1,size(intersections) + dsps = dsp_line_from_dsp(intersections(i), pt_hpl, pt_bub) + end do + end block dsp_line_hpl_bub call pt_three_phase_from_intersection(& pt_hpl, pt_bub, intersections, & pt_dew_3, pt_bub_3 & From b36fbab5bef537b2fd530a492d0192f1f1a4c7e6 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:27:26 -0300 Subject: [PATCH 28/29] calculation of self intersection on Px lines, should add bubble --- app/main.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/app/main.f90 b/app/main.f90 index eef673b..ef168db 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -278,7 +278,7 @@ subroutine px_envelopes ! Look for crossings ! ------------------------------------------------------------------------ inter = intersection(px_dew%alpha, px_dew%p, px_bub%alpha, px_bub%p) - self_inter = intersection(px_dew%alpha, px_dew%p) + self_inter = intersection(px_bub%alpha, px_bub%p) print *, style_bold // "Px Intersections: " // style_reset, size(inter) print *, style_bold // "Px Self-Intersections: " // style_reset, size(self_inter) ! ======================================================================== @@ -307,7 +307,7 @@ subroutine px_envelopes !TODO: Add a check if one of the previous lines already ! in this DSP px_branch_3 = px_three_phase_from_inter(& - self_inter(i), px_dew, px_dew & + self_inter(i), px_bub, px_bub & ) end do end if From b9bcf2f88a7d085d9c20361eedf22c87851f63fc Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Wed, 15 Nov 2023 10:27:35 -0300 Subject: [PATCH 29/29] tolerances for PxfromPT --- app/main.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/app/main.f90 b/app/main.f90 index ef168db..b563c1b 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -268,7 +268,7 @@ subroutine px_envelopes px_bub = px_two_phase_from_pt(t_inj, pt_bub, t_tol=5.0_pr) print *, blue // "Running Dew" // style_reset - px_dew = px_two_phase_from_pt(t_inj, pt_dew, t_tol=15.0_pr) + px_dew = px_two_phase_from_pt(t_inj, pt_dew, t_tol=5.0_pr) print *, blue // "Running HPLL" // style_reset px_hpl = px_hpl_line(0.99_pr, maxval(px_bub%p))