DSP lines

This is a big pull request. Adding the calculation of new lines (DSP lines). Which start from a single DSP obtained from PT lines

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
@@ -106,18 +109,18 @@ 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 &
-         )
+      )
       ! ========================================================================
 
       ! ========================================================================
       !  Dew/AOP envelopes
       ! ------------------------------------------------------------------------
-      t = 200
+      t = 300
       p = p_wilson(z, t)
       do while (p > 0.1)
          t = t - 5
@@ -131,19 +134,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
       ! ========================================================================
 
       ! ========================================================================
@@ -176,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 &
@@ -193,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 &
@@ -238,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))
@@ -248,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)
       ! ========================================================================
@@ -277,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

fpm.toml

@@ -6,7 +6,7 @@ maintainer = "[email protected]"
 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" }

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), &
@@ -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