From c81349c6aa2817c59d61a3f9222835dc7802108a Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 11:05:18 -0300 Subject: [PATCH 01/15] examples --- example/tutorials/basics.f90 | 25 +++++++++++---------- example/tutorials/new_alpha_function.f90 | 2 -- example/tutorials/phase_envelope.f90 | 28 +++++++++++++++++++++--- example/tutorials/saturation_points.f90 | 10 ++++----- 4 files changed, 44 insertions(+), 21 deletions(-) diff --git a/example/tutorials/basics.f90 b/example/tutorials/basics.f90 index 0b6e2b653..70554bcf5 100644 --- a/example/tutorials/basics.f90 +++ b/example/tutorials/basics.f90 @@ -1,18 +1,21 @@ program basics - use yaeos ! First the module is imported (used) + use yaeos ! First the module is imported (used). This will bring all the + ! types and functions defined in the module to this program. implicit none ! =========================================================================== - ! In yaeos all residual Helmholtz models are based on the structure - ! `ArModel` and are derivatives from it. allocatable means that - ! this variable can be redefined as any other kind of variable that derives - ! from `ArModel`, for example `CubicEoS`. + + ! In yaeos, all residual Helmholtz models are based on the structure + ! `ArModel` and are derivatives from it. + ! Defining the model variable as an allocatable ArModel means that this + ! variable can be redefined as any other kind of variable that derives from + ! `ArModel`, for example `CubicEoS`. class(ArModel), allocatable :: model - ! Number of components - integer, parameter :: nc=2 + integer, parameter :: nc=2 !! In this case we define a constant number of + !! components. ! Used input variables - real(pr) :: n(nc), tc(nc), pc(nc), w(nc) + real(pr) :: n(nc), Tc(nc), Pc(nc), w(nc) ! Output variables real(pr) :: P, dPdV @@ -32,10 +35,10 @@ program basics n = [0.3, 0.7] ! Pressure calculation - call model%pressure(n, v=2.5_pr, T=150._pr, P=P) + call model%pressure(n, V=2.5_pr, T=150._pr, P=P) print *, "P: ", P ! Derivatives can also be calculated when included as optional arguments! - call model%pressure(n, v=2.5_pr, T=150._pr, P=P, dPdV=dPdV) + call model%pressure(n, V=2.5_pr, T=150._pr, P=P, dPdV=dPdV) print *, "dPdV: ", dPdV -end program basics +end program basics \ No newline at end of file diff --git a/example/tutorials/new_alpha_function.f90 b/example/tutorials/new_alpha_function.f90 index c222936c9..ec261f007 100644 --- a/example/tutorials/new_alpha_function.f90 +++ b/example/tutorials/new_alpha_function.f90 @@ -75,6 +75,4 @@ program new_alpha_example call eos%pressure(n, V, T, P=P) print *, "Peng-Robinson76-MC:", P - - end program new_alpha_example diff --git a/example/tutorials/phase_envelope.f90 b/example/tutorials/phase_envelope.f90 index a2d82b7b5..381eb2aaa 100644 --- a/example/tutorials/phase_envelope.f90 +++ b/example/tutorials/phase_envelope.f90 @@ -3,8 +3,8 @@ program phase_envelope use yaeos__example_tools, only: methane_butane_pr76 implicit none - class(ArModel), allocatable :: model ! model - type(EquilibriumState) :: bubble ! + class(ArModel), allocatable :: model ! Model + type(EquilibriumState) :: bubble ! Saturation point type(PTEnvel2) :: envelope integer, parameter :: nc=2 @@ -18,7 +18,7 @@ program phase_envelope ! Calculate a bubble point at 100K to initialize the rest of the phase ! envelope calculation - bubble = saturation_pressure(model, z, T=100._pr, kind="bubble", p0=40._pr) + bubble = saturation_pressure(model, z, T=150._pr, kind="bubble") ! Calculate the whole phase envelope using the information from the converged ! dew point @@ -26,4 +26,26 @@ program phase_envelope ! Write the resulting phase envelope to the screen write(*, *) envelope + + ! Write the resulting phase envelope to a file + open(1, file="phase_envelope.dat") + write(1, *) envelope + close(1) + + ! It is also possible to initialize the phase_envelope with a manually + ! defined point + manual: block + real(pr) :: k(nc), T, P + T = 150 + P = 1 + k = k_wilson(model, T=T, P=P) + + ! Define manually an initial point + bubble%kind = "bubble" + bubble%x = z + bubble%y = k*z + bubble%beta = 0 + + envelope = pt_envelope_2ph(model, z, bubble) + end block manual end program phase_envelope \ No newline at end of file diff --git a/example/tutorials/saturation_points.f90 b/example/tutorials/saturation_points.f90 index 26b8c1fe8..e099fd74c 100644 --- a/example/tutorials/saturation_points.f90 +++ b/example/tutorials/saturation_points.f90 @@ -24,7 +24,7 @@ program saturation write(*, *) "Bubble pressure:" T = 150 sat_point = saturation_pressure(model, n, T=T, kind="bubble") - write (*, *) "kind, T, P: ", sat_point + write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P write (*, *) "x: ", sat_point%x write (*, *) "y: ", sat_point%y @@ -34,7 +34,7 @@ program saturation write(*, *) "" write (*, *) "Bubble temperature:" sat_point = saturation_temperature(model, n, P=15._pr, kind="bubble") - write (*, *) "kind, T, P: ", sat_point + write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P write (*, *) "x: ", sat_point%x write (*, *) "y: ", sat_point%y @@ -44,7 +44,7 @@ program saturation write(*, *) "" write(*, *) "Dew pressure:" sat_point = saturation_pressure(model, n, T=150._pr, kind="dew") - write (*, *) "kind, T, P: ", sat_point + write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P write (*, *) "x: ", sat_point%x write (*, *) "y: ", sat_point%y @@ -53,8 +53,8 @@ program saturation ! -------------------------------------------------------------------------- write(*, *) "" write (*, *) "Dew temperature:" - sat_point = saturation_temperature(model, n, P=15._pr, kind="dew") - write (*, *) "kind, T, P: ", sat_point + sat_point = saturation_temperature(model, n, P=15._pr, kind="dew", t0=330._pr) + write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P write (*, *) "x: ", sat_point%x write (*, *) "y: ", sat_point%y end program \ No newline at end of file From 7fb1521b55d211abf5ece1c891daf231ce7213bf Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 11:06:00 -0300 Subject: [PATCH 02/15] expose k_wilson correlations --- src/phase_equilibria/phase_equilibria.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/phase_equilibria/phase_equilibria.f90 b/src/phase_equilibria/phase_equilibria.f90 index 29d1c8bdd..c8a415ece 100644 --- a/src/phase_equilibria/phase_equilibria.f90 +++ b/src/phase_equilibria/phase_equilibria.f90 @@ -5,5 +5,6 @@ module yaeos__equilibria saturation_pressure, saturation_temperature use yaeos__phase_equilibria_boundaries_phase_envelopes_pt, only:& CriticalPoint, PTEnvel2, pt_envelope_2ph + use yaeos__phase_equilibria_auxiliar, only: k_wilson implicit none end module \ No newline at end of file From b69842d0ec5ca106d74cc01f2c3c712e139861fe Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 12:25:08 -0300 Subject: [PATCH 03/15] thermoprops example --- example/tutorials/basics.f90 | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/example/tutorials/basics.f90 b/example/tutorials/basics.f90 index 70554bcf5..f73febd00 100644 --- a/example/tutorials/basics.f90 +++ b/example/tutorials/basics.f90 @@ -41,4 +41,30 @@ program basics ! Derivatives can also be calculated when included as optional arguments! call model%pressure(n, V=2.5_pr, T=150._pr, P=P, dPdV=dPdV) print *, "dPdV: ", dPdV + + thermoproperties: block + real(pr) :: Cpr, Cvr + real(pr) :: Gr, GrT, GrV, Grn(nc) + real(pr) :: Hr, HrT, HrV, Hrn(nc) + real(pr) :: lnPhi(nc), dlnPhidT(nc), dlnPhidP(nc), dlnPhidn(nc, nc) + + real(pr) :: n(nc), P, V, T + + n = [0.3_pr, 0.7_pr] + V = 3.5_pr + P = 5.0_pr + call model%Cp_residual_vt(n, V=2.5_pr, T=150._pr, Cp=Cpr) + call model%Cv_residual_vt(n, V=2.5_pr, T=150._pr, Cv=Cvr) + call model%gibbs_residual_vt(n, V=2.5_pr, T=150._pr, Gr=Gr, GrT=GrT, GrV=GrV, Grn=Grn) + call model%enthalpy_residual_vt(n, V=2.5_pr, T=150._pr, Hr=Hr, HrT=HrT, HrV=HrV, Hrn=Hrn) + call model%lnphi_pt(& + n, P=P, T=T, root_type="stable", & + lnPhi=lnPhi, dlnPhidT=dlnphidT, dlnPhidP=dlnPhidP, dlnPhidn=dlnPhidn & + ) + call model%lnphi_vt(& + n, V=V, T=T, & + lnPhi=lnPhi, dlnPhidT=dlnphidT, dlnPhidP=dlnPhidP, dlnPhidn=dlnPhidn & + ) + call model%volume(n, P, T, root_type="vapor", V=V) + end block thermoproperties end program basics \ No newline at end of file From cc9cf5baa3404d3b87dfa00412e298018adfbd97 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 12:26:06 -0300 Subject: [PATCH 04/15] doc --- example/tutorials/basics.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/example/tutorials/basics.f90 b/example/tutorials/basics.f90 index f73febd00..ded5886aa 100644 --- a/example/tutorials/basics.f90 +++ b/example/tutorials/basics.f90 @@ -43,6 +43,7 @@ program basics print *, "dPdV: ", dPdV thermoproperties: block + !! Al the bulk thermodynamic properties currently available real(pr) :: Cpr, Cvr real(pr) :: Gr, GrT, GrV, Grn(nc) real(pr) :: Hr, HrT, HrV, Hrn(nc) From 1e36844dca21c21939f80e8ca6e66c49bf8cfba9 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 12:35:51 -0300 Subject: [PATCH 05/15] examples names --- example/tutorials/{basics.f90 => 1_basics.f90} | 2 +- .../{saturation_points.f90 => 2_saturation_points.f90} | 2 ++ fpm.toml | 4 ++-- src/models/residual_helmholtz/cubic/generic_cubic.f90 | 2 +- 4 files changed, 6 insertions(+), 4 deletions(-) rename example/tutorials/{basics.f90 => 1_basics.f90} (97%) rename example/tutorials/{saturation_points.f90 => 2_saturation_points.f90} (94%) diff --git a/example/tutorials/basics.f90 b/example/tutorials/1_basics.f90 similarity index 97% rename from example/tutorials/basics.f90 rename to example/tutorials/1_basics.f90 index ded5886aa..a3fb7653c 100644 --- a/example/tutorials/basics.f90 +++ b/example/tutorials/1_basics.f90 @@ -66,6 +66,6 @@ program basics n, V=V, T=T, & lnPhi=lnPhi, dlnPhidT=dlnphidT, dlnPhidP=dlnPhidP, dlnPhidn=dlnPhidn & ) - call model%volume(n, P, T, root_type="vapor", V=V) + call model%volume(n, P=P, T=T, root_type="vapor", V=V) end block thermoproperties end program basics \ No newline at end of file diff --git a/example/tutorials/saturation_points.f90 b/example/tutorials/2_saturation_points.f90 similarity index 94% rename from example/tutorials/saturation_points.f90 rename to example/tutorials/2_saturation_points.f90 index e099fd74c..391c4fa03 100644 --- a/example/tutorials/saturation_points.f90 +++ b/example/tutorials/2_saturation_points.f90 @@ -14,6 +14,8 @@ program saturation real(pr) :: n(2), T + ! This function returns a predefined model to simplify the step. How it + ! works can be seen in the `example_tools.f90` file. model = methane_butane_pr76() n = [2.5, 6.7] diff --git a/fpm.toml b/fpm.toml index d1906b24b..080bf4291 100644 --- a/fpm.toml +++ b/fpm.toml @@ -39,12 +39,12 @@ main = "demo.f90" [[example]] name = "basics" source-dir = "example/tutorials" -main = "basics.f90" +main = "1_basics.f90" [[example]] name = "saturation_points" source-dir = "example/tutorials" -main = "saturation_points.f90" +main = "2_saturation_points.f90" [[example]] name = "phase_envelope" diff --git a/src/models/residual_helmholtz/cubic/generic_cubic.f90 b/src/models/residual_helmholtz/cubic/generic_cubic.f90 index ea80156b6..6e9812612 100644 --- a/src/models/residual_helmholtz/cubic/generic_cubic.f90 +++ b/src/models/residual_helmholtz/cubic/generic_cubic.f90 @@ -345,7 +345,7 @@ subroutine volume(eos, n, P, T, V, root_type) real(pr) :: a(size(n)), dadt(size(n)), dadt2(size(n)) real(pr) :: totn - call volume_michelsen(eos, n, P, T, V, root_type) + call volume_michelsen(eos, n=n, P=P, T=T, V=V, root_type=root_type) return totn = sum(n) From 10e425462194f90262b22e2c122ae0e6d360acd6 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 12:35:58 -0300 Subject: [PATCH 06/15] quadsendo --- src/math/linalg.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/math/linalg.f90 b/src/math/linalg.f90 index 82523eb85..897f3022a 100644 --- a/src/math/linalg.f90 +++ b/src/math/linalg.f90 @@ -107,22 +107,22 @@ subroutine cubic_roots_rosendo(parameters, real_roots, complex_roots, flag) complex(pr), intent(out) :: complex_roots(3) integer, intent(out) :: flag - real(pr) :: d1, d2, d3, Q, R, A, B, theta, alp, bet, gam + real(16) :: d1, d2, d3, Q, R, A, B, theta, alp, bet, gam integer :: i d1 = parameters(2) / parameters(1) d2 = parameters(3) / parameters(1) d3 = parameters(4) / parameters(1) - Q = (d1**2 - 3*d2) / 9.0_pr - R = (2*d1**3 - 9*d1*d2 + 27*d3) / 54.0_pr + Q = (d1**2 - 3*d2) / 9.0_16 + R = (2*d1**3 - 9*d1*d2 + 27*d3) / 54.0_16 if (R**2 <= Q**3) then theta = acos(R / sqrt(Q**3)) - real_roots(1) = -2 * sqrt(Q) * cos(theta / 3.0_pr) - d1 / 3.0_pr - real_roots(2) = -2 * sqrt(Q) * cos((theta + 2 * pi) / 3.0_pr) - d1 / 3.0_pr - real_roots(3) = -2 * sqrt(Q) * cos((theta - 2 * pi) / 3.0_pr) - d1 / 3.0_pr + real_roots(1) = -2 * sqrt(Q) * cos(theta / 3.0_16) - d1 / 3.0_16 + real_roots(2) = -2 * sqrt(Q) * cos((theta + 2 * pi) / 3.0_16) - d1 / 3.0_16 + real_roots(3) = -2 * sqrt(Q) * cos((theta - 2 * pi) / 3.0_16) - d1 / 3.0_16 ! Correction?? ! do i=1,100 @@ -134,16 +134,16 @@ subroutine cubic_roots_rosendo(parameters, real_roots, complex_roots, flag) call sort(real_roots) flag = -1 else - A = - sign((abs(R) + sqrt(R**2 - Q**3))**(1.0_pr/3.0_pr), R) + A = - sign((abs(R) + sqrt(R**2 - Q**3))**(1.0_16/3.0_16), R) if (abs(A) < 1e-6) then - A = 0.0_pr - B = 0.0_pr + A = 0.0_16 + B = 0.0_16 else B = Q / A end if - real_roots = (A + B) - d1 / 3.0_pr + real_roots = (A + B) - d1 / 3.0_16 flag = 1 end if end subroutine From 1d61856c9572038e8037fefa380c2ea29d9b726e Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 13:44:35 -0300 Subject: [PATCH 07/15] flash example --- example/tutorials/3_phase_split.f90 | 46 +++++++++++++++++++++ fpm.toml | 5 +++ src/phase_equilibria/saturations_points.f90 | 13 ++++++ 3 files changed, 64 insertions(+) create mode 100644 example/tutorials/3_phase_split.f90 diff --git a/example/tutorials/3_phase_split.f90 b/example/tutorials/3_phase_split.f90 new file mode 100644 index 000000000..700afbd06 --- /dev/null +++ b/example/tutorials/3_phase_split.f90 @@ -0,0 +1,46 @@ +program flasher + !! Example program to make flash calculations with a PengRobinson76 + !! EoS model. + + ! Import the relevant + use yaeos, only: pr, EquilibriumState, flash, PengRobinson76, ArModel, k_wilson + implicit none + + ! Variables definition: + class(ArModel), allocatable :: model !! Model to use + type(EquilibriumState) :: flash_result !! Result of Flash calculation + + real(pr) :: tc(2), pc(2), w(2) + real(pr) :: n(2), t, p, k0(2) + integer :: iter + + ! Code starts here: + print *, "FLASH EXAMPLE:" + print *, "Methane/Butane mixture" + + ! Methane/ Butane mixture + n = [0.4, 0.6] ! Composition + tc = [190.564, 425.12] ! Critical temperatures + pc = [45.99, 37.96] ! Critical pressures + w = [0.0115478, 0.200164] ! Acentric factors + + ! Use the PengRobinson76 model + model = PengRobinson76(tc, pc, w) + + ! Set pressure and temperatures + P = 60 + T = 294 + + ! Calculate flashes + flash_result = flash(model, n, t=t, p_spec=p, iters=iter) + write(*, *) flash_result + + ! K-wilson factors for initialization + k0 = k_wilson(model, T=T, P=P) + flash_result = flash(model, n, t=t, p_spec=p, k0=k0, iters=iter) + write(*, *) flash_result + + ! Specify volume + flash_result = flash(model, n, t=t, V_spec=1.0_pr, iters=iter) + write(*, *) flash_result +end program diff --git a/fpm.toml b/fpm.toml index 080bf4291..38917f84a 100644 --- a/fpm.toml +++ b/fpm.toml @@ -46,6 +46,11 @@ name = "saturation_points" source-dir = "example/tutorials" main = "2_saturation_points.f90" +[[example]] +name = "phase_split" +source-dir = "example/tutorials" +main = "3_phase_split.f90" + [[example]] name = "phase_envelope" source-dir = "example/tutorials" diff --git a/src/phase_equilibria/saturations_points.f90 b/src/phase_equilibria/saturations_points.f90 index ef8377a64..4f0321fe4 100644 --- a/src/phase_equilibria/saturations_points.f90 +++ b/src/phase_equilibria/saturations_points.f90 @@ -150,9 +150,12 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0 real(pr) :: f, step integer :: its, iterations + logical :: is_incipient(size(n)) + ! ======================================================================= ! Handle arguments ! ----------------------------------------------------------------------- + is_incipient = .true. z = n/sum(n) if (present (t0)) then t = t0 @@ -185,6 +188,11 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0 where (z == 0) k = 0 end where + + where (y == 0) + is_incipient = .false. + end where + ! ======================================================================== ! ======================================================================== @@ -192,6 +200,10 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0 ! ------------------------------------------------------------------------ do its=1, iterations y = k*z + where (.not. is_incipient) + y = 0 + endwhere + call model%lnphi_pt(y, P, T, vy, incipient, lnPhi=lnfug_y, dlnphidt=dlnphi_dt_y) call model%lnphi_pt(z, P, T, vz, main, lnPhi=lnfug_z, dlnphidt=dlnphi_dt_z) @@ -204,6 +216,7 @@ type(EquilibriumState) function saturation_temperature(model, n, p, kind, t0, y0 end do t = t - step + if (abs(step) < tol .and. abs(f) < tol) exit end do ! ======================================================================== From 16dae25157fccf89be862a8c3ab47d8dd228e192 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 13:50:43 -0300 Subject: [PATCH 08/15] renamed phase envel example --- example/tutorials/{phase_envelope.f90 => 4_phase_envelope.f90} | 0 fpm.toml | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename example/tutorials/{phase_envelope.f90 => 4_phase_envelope.f90} (100%) diff --git a/example/tutorials/phase_envelope.f90 b/example/tutorials/4_phase_envelope.f90 similarity index 100% rename from example/tutorials/phase_envelope.f90 rename to example/tutorials/4_phase_envelope.f90 diff --git a/fpm.toml b/fpm.toml index 38917f84a..63bad66d4 100644 --- a/fpm.toml +++ b/fpm.toml @@ -54,7 +54,7 @@ main = "3_phase_split.f90" [[example]] name = "phase_envelope" source-dir = "example/tutorials" -main = "phase_envelope.f90" +main = "4_phase_envelope.f90" [[example]] name = "new_alpha_function" From 16e15ad1d4993e18394c11b545efaae8257e42f5 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 14:03:50 -0300 Subject: [PATCH 09/15] moved to another folder --- example/{tutorials => basics}/1_basics.f90 | 0 .../2_saturation_points.f90 | 0 .../{tutorials => basics}/3_phase_split.f90 | 0 .../4_phase_envelope.f90 | 0 .../pure_psat.f90 => basics/5_pure_psat.f90} | 20 +- example/{tutorials => basics}/cubic_eos.f90 | 0 .../{tutorials => basics}/example_tools.f90 | 0 example/{tutorials => basics}/huron_vidal.f90 | 0 .../new_alpha_function.f90 | 0 example/extra/benchmark.f90 | 4 +- example/extra/demo.f90 | 24 +- example/extra/tape_pr76.f90 | 49 - example/extra/taperobinson.f90 | 1354 ----------------- fpm.toml | 18 +- 14 files changed, 35 insertions(+), 1434 deletions(-) rename example/{tutorials => basics}/1_basics.f90 (100%) rename example/{tutorials => basics}/2_saturation_points.f90 (100%) rename example/{tutorials => basics}/3_phase_split.f90 (100%) rename example/{tutorials => basics}/4_phase_envelope.f90 (100%) rename example/{extra/pure_psat.f90 => basics/5_pure_psat.f90} (81%) rename example/{tutorials => basics}/cubic_eos.f90 (100%) rename example/{tutorials => basics}/example_tools.f90 (100%) rename example/{tutorials => basics}/huron_vidal.f90 (100%) rename example/{tutorials => basics}/new_alpha_function.f90 (100%) delete mode 100644 example/extra/tape_pr76.f90 delete mode 100644 example/extra/taperobinson.f90 diff --git a/example/tutorials/1_basics.f90 b/example/basics/1_basics.f90 similarity index 100% rename from example/tutorials/1_basics.f90 rename to example/basics/1_basics.f90 diff --git a/example/tutorials/2_saturation_points.f90 b/example/basics/2_saturation_points.f90 similarity index 100% rename from example/tutorials/2_saturation_points.f90 rename to example/basics/2_saturation_points.f90 diff --git a/example/tutorials/3_phase_split.f90 b/example/basics/3_phase_split.f90 similarity index 100% rename from example/tutorials/3_phase_split.f90 rename to example/basics/3_phase_split.f90 diff --git a/example/tutorials/4_phase_envelope.f90 b/example/basics/4_phase_envelope.f90 similarity index 100% rename from example/tutorials/4_phase_envelope.f90 rename to example/basics/4_phase_envelope.f90 diff --git a/example/extra/pure_psat.f90 b/example/basics/5_pure_psat.f90 similarity index 81% rename from example/extra/pure_psat.f90 rename to example/basics/5_pure_psat.f90 index 90ed05307..c0cef2028 100644 --- a/example/extra/pure_psat.f90 +++ b/example/basics/5_pure_psat.f90 @@ -1,16 +1,17 @@ !! Program to calculate the vapor pressure of pure components -!! - module pure_psat !! Module used to calculate the saturation pressure of pure components at !! a given temperature. use yaeos contains real(pr) function Psat(eos, ncomp, T) - use yaeos__math, only: solve_system - class(ArModel), intent(in) :: eos - integer, intent(in) :: ncomp - real(pr), intent(in) :: T + !! Calculation of saturation pressure of a pure component using the + !! secant method. + class(ArModel), intent(in) :: eos !! Model that will be used + integer, intent(in) :: ncomp + !! Number of component in the mixture from which the saturation pressure + !! will be calculated + real(pr), intent(in) :: T !! Temperature [K] real(pr) :: P1, P2 real(pr) :: f1, f2 @@ -30,7 +31,6 @@ real(pr) function Psat(eos, ncomp, T) P1 = P2 P2 = Psat end do - contains real(pr) function diff(P) real(pr), intent(in) :: P @@ -47,11 +47,13 @@ program main use yaeos use forsus, only: Substance, forsus_default_dir, forsus_dir use pure_psat, only: Psat + implicit none + integer, parameter :: nc=2 type(CubicEoS) :: eos type(Substance) :: sus(nc) - real(pr) :: n(nc), T, f + real(pr) :: n(nc), T integer :: i, j forsus_dir = "build/dependencies/forsus/" // forsus_default_dir @@ -62,6 +64,8 @@ program main sus%critical%critical_pressure%value/1e5,& sus%critical%acentric_factor%value & ) + + T = 273.15_pr + 50 do i=273+90, nint(maxval(sus%critical%critical_temperature%value)) T = real(i,pr) diff --git a/example/tutorials/cubic_eos.f90 b/example/basics/cubic_eos.f90 similarity index 100% rename from example/tutorials/cubic_eos.f90 rename to example/basics/cubic_eos.f90 diff --git a/example/tutorials/example_tools.f90 b/example/basics/example_tools.f90 similarity index 100% rename from example/tutorials/example_tools.f90 rename to example/basics/example_tools.f90 diff --git a/example/tutorials/huron_vidal.f90 b/example/basics/huron_vidal.f90 similarity index 100% rename from example/tutorials/huron_vidal.f90 rename to example/basics/huron_vidal.f90 diff --git a/example/tutorials/new_alpha_function.f90 b/example/basics/new_alpha_function.f90 similarity index 100% rename from example/tutorials/new_alpha_function.f90 rename to example/basics/new_alpha_function.f90 diff --git a/example/extra/benchmark.f90 b/example/extra/benchmark.f90 index 6cd8661f0..db1f227a9 100644 --- a/example/extra/benchmark.f90 +++ b/example/extra/benchmark.f90 @@ -2,7 +2,7 @@ module bench use yaeos, only: pr, R, Substances, AlphaSoave, CubicEoS, & QMR, PengRobinson76, ArModel use hyperdual_pr76, only: PR76, setup_adiff_pr76 => setup - use TapeRobinson, only: setup_taperobinson => setup_model + ! use TapeRobinson, only: setup_taperobinson => setup_model implicit none real(pr), allocatable :: z(:), tc(:), pc(:), w(:), kij(:,:), lij(:,:) @@ -40,7 +40,7 @@ subroutine yaeos__run(n, dn, f_p, model_name) case ("Adiff PR76") model = setup_adiff_pr76(tc, pc, w, kij, lij) case ("Tape PR76") - model = setup_taperobinson(tc, pc, w, kij, lij) + ! model = setup_taperobinson(tc, pc, w, kij, lij) end select v = 1.0_pr diff --git a/example/extra/demo.f90 b/example/extra/demo.f90 index 0ba70a481..7d256fae0 100644 --- a/example/extra/demo.f90 +++ b/example/extra/demo.f90 @@ -1,18 +1,18 @@ program examples - use bench, only: benchmarks => main + ! use bench, only: benchmarks => main use hyperdual_pr76, only: run_hyperdual_pr76 => main use flashing, only: run_flashes => main - use TapeRobinson, only: run_tape_pr76 => main use tape_nrtl, only: run_tape_nrtl => main - ! print *, "Running Tapenade generated NRTL model" - ! call run_tape_nrtl - ! print *, "Running Tapenade generated PR76" - ! call run_tape_pr76 - ! print *, "Running Hyperdual generated PR76" - ! call run_hyperdual_pr76 - print *, "Running bencharks O(f(N))" - call benchmarks - ! print *, "Flash example" - ! call run_flashes + print *, "=========================================" + print *, "YAEOS DEMO" + print *, "Running Tapenade generated NRTL model" + call run_tape_nrtl + print *, "Running Hyperdual generated PR76" + call run_hyperdual_pr76 + ! print *, "Running bencharks O(f(N))" + ! call benchmarks + print *, "Flash example" + call run_flashes + print *, "=========================================" end program \ No newline at end of file diff --git a/example/extra/tape_pr76.f90 b/example/extra/tape_pr76.f90 deleted file mode 100644 index 88d9c9649..000000000 --- a/example/extra/tape_pr76.f90 +++ /dev/null @@ -1,49 +0,0 @@ -MODULE TapeRobinson - use tapenade_pr, only: SETUP_MODEL -contains - subroutine main() - use yaeos, only: ArModel - use yaeos__constants, only: pr - class(ArModel), allocatable :: model - integer, parameter :: n=2 - real(pr) :: z(n), tc(n), pc(n), w(n), kij(n,n), lij(n,n) - - real(pr) :: v, t - - real(pr) :: ar, arv, arv2, art, art2, artv - real(pr) :: arn(n), arvn(n), artn(n), arn2(n,n) - - z = [0.3_pr, 0.7_pr] - tc = [190._pr, 310._pr] - pc = [14._pr, 30._pr] - w = [0.001_pr, 0.03_pr] - - kij = reshape([0._pr, 0.1_pr, 0.1_pr, 0._pr], [n,n]) - lij = kij / 2._pr - - v = 1 - t = 150 - - model = SETUP_MODEL(tc, pc, w, kij, lij) - call model%residual_helmholtz(& - z, V, T, Ar=Ar, ArV=ArV, ArV2=ArV2, ArT=ArT, ArTV=ArTV, & - ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2 & - ) - print *, "Ar: ", ar - - print *, "ArV: ", arV - print *, "ArT: ", arT - - print *, "ArT2: ", arT2 - print *, "ArV2: ", ArV2 - - print *, "ArTV: ", ArTV - - print *, "Arn: ", Arn - - print *, "ArVn: ", ArVn - print *, "ArTn: ", ArTn - - print *, "Arn2: ", Arn2 - end subroutine main -end module TapeRobinson diff --git a/example/extra/taperobinson.f90 b/example/extra/taperobinson.f90 deleted file mode 100644 index 72dba51ad..000000000 --- a/example/extra/taperobinson.f90 +++ /dev/null @@ -1,1354 +0,0 @@ -! Generated by TAPENADE (INRIA, Ecuador team) -! Tapenade 3.16 (develop) - 13 Sep 2023 12:36 -! -! Generated by TAPENADE (INRIA, Ecuador team) -! Tapenade 3.16 (develop) - 13 Sep 2023 12:36 -! -! Generated by TAPENADE (INRIA, Ecuador team) -! Tapenade 3.16 (develop) - 13 Sep 2023 12:36 -! -! Generated by TAPENADE (INRIA, Ecuador team) -! Tapenade 3.16 (develop) - 13 Sep 2023 12:36 -! -! Generated by TAPENADE (INRIA, Ecuador team) -! Tapenade 3.16 (develop) - 13 Sep 2023 12:36 -! -MODULE TAPENADE_PR - USE YAEOS__CONSTANTS, ONLY : pr, r - USE YAEOS__TAPENADE_AR_API, ONLY : armodeltapenade - use yaeos__tapenade_interfaces - IMPLICIT NONE -type, extends(ArModelTapenade) :: TPR76 - REAL(pr), ALLOCATABLE :: kij(:, :), lij(:, :) - REAL(pr), ALLOCATABLE :: ac(:), b(:), k(:) - REAL(pr), ALLOCATABLE :: tc(:), pc(:), w(:) - REAL(pr) :: r=0.08314472_pr - REAL(pr) :: del1=1.+SQRT(2.) - REAL(pr) :: del2=1.-SQRT(2.) - contains - procedure :: ar - procedure :: ar_d - procedure :: ar_b - procedure :: ar_d_b - procedure :: ar_d_d - procedure :: v0 => VOLUME_INITALIZER - END TYPE TPR76 - -CONTAINS - TYPE(TPR76) FUNCTION SETUP_MODEL(tc, pc, w, kij, lij) - IMPLICIT NONE - REAL(pr), INTENT(IN) :: tc(:) - REAL(pr), INTENT(IN) :: pc(:) - REAL(pr), INTENT(IN) :: w(:) - REAL(pr), INTENT(IN) :: kij(:, :) - REAL(pr), INTENT(IN) :: lij(:, :) - - setup_model%tc = tc - setup_model%pc = pc - setup_model%w = w - setup_model%ac = 0.45723553*r**2*tc**2/pc - setup_model%b = 0.07779607*r*tc/pc - setup_model%k = 0.37464 + 1.54226*w - 0.26993*w**2 - setup_model%kij = kij - setup_model%lij = lij - END FUNCTION SETUP_MODEL - -! Differentiation of ar_d_d in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval arvaldd arvald arvald0 -! with respect to varying inputs: t v -! RW status of diff variables: t:in v:in arval:out arvaldd:out -! arvald:out arvald0:out -! Differentiation of ar_d in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval arvald -! with respect to varying inputs: t v -! RW status of diff variables: t:in v:in arval:out arvald:out -! Differentiation of ar in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval -! with respect to varying inputs: n t v -! RW status of diff variables: n:in t:in v:in arval:out - SUBROUTINE AR_D_D_D(model, n, nd, v, vd1, vd0, vd, t, td1, td0, td, & -& arval, arvald1, arvald0, arvald0d, arvald, arvaldd0, arvaldd, & -& arvalddd) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:), v, t - REAL(pr), INTENT(IN) :: vd1, td1 - REAL(pr), INTENT(IN) :: vd0, td0 - REAL(pr), INTENT(IN) :: nd(:), vd, td - REAL(pr), INTENT(OUT) :: arval - REAL(pr), INTENT(OUT) :: arvald1 - REAL(pr), INTENT(OUT) :: arvald0 - REAL(pr), INTENT(OUT) :: arvald0d - REAL(pr), INTENT(OUT) :: arvald - REAL(pr), INTENT(OUT) :: arvaldd0 - REAL(pr), INTENT(OUT) :: arvaldd - REAL(pr), INTENT(OUT) :: arvalddd - REAL(pr) :: amix, a(SIZE(n)), ai(SIZE(n)), z2(SIZE(n)), nij - REAL(pr) :: amixd1, ad1(SIZE(n)) - REAL(pr) :: amixd0, ad0(SIZE(n)) - REAL(pr) :: amixd0d, ad0d(SIZE(n)) - REAL(pr) :: amixd, ad(SIZE(n)), nijd - REAL(pr) :: amixdd0, add0(SIZE(n)) - REAL(pr) :: amixdd, add(SIZE(n)) - REAL(pr) :: amixddd, addd(SIZE(n)) - REAL(pr) :: bmix - REAL(pr) :: bmixd - REAL(pr) :: b_v - REAL(pr) :: b_vd1 - REAL(pr) :: b_vd0 - REAL(pr) :: b_vd0d - REAL(pr) :: b_vd - REAL(pr) :: b_vdd0 - REAL(pr) :: b_vdd - REAL(pr) :: b_vddd - REAL(pr) :: aij(SIZE(n), SIZE(n)), bij(SIZE(n), SIZE(n)) - INTEGER :: i, j - REAL(pr) :: tc(SIZE(n)), pc(SIZE(n)), w(SIZE(n)) - REAL(pr) :: ac(SIZE(n)), b(SIZE(n)), k(SIZE(n)), r - REAL(pr) :: del1, del2 - REAL(pr) :: kij(SIZE(n), SIZE(n)), lij(SIZE(n), SIZE(n)) - INTRINSIC SQRT - INTRINSIC SUM - INTRINSIC LOG - INTRINSIC SIZE - REAL(pr), DIMENSION(SIZE(n)) :: arg1 - REAL(pr), DIMENSION(SIZE(n)) :: arg1d1 - REAL(pr), DIMENSION(SIZE(n)) :: arg1d0 - REAL(pr), DIMENSION(SIZE(n)) :: arg1d0d - REAL(pr), DIMENSION(SIZE(n)) :: arg1d - REAL(pr), DIMENSION(SIZE(n)) :: arg1dd0 - REAL(pr), DIMENSION(SIZE(n)) :: arg1dd - REAL(pr), DIMENSION(SIZE(n)) :: arg1ddd - REAL(pr), DIMENSION(SIZE(n)) :: result1 - REAL(pr), DIMENSION(SIZE(n)) :: result1d1 - REAL(pr), DIMENSION(SIZE(n)) :: result1d0 - REAL(pr), DIMENSION(SIZE(n)) :: result1d0d - REAL(pr), DIMENSION(SIZE(n)) :: result1d - REAL(pr), DIMENSION(SIZE(n)) :: result1dd0 - REAL(pr), DIMENSION(SIZE(n)) :: result1dd - REAL(pr), DIMENSION(SIZE(n)) :: result1ddd - REAL(pr), DIMENSION(SIZE(n)) :: arg2 - REAL(pr), DIMENSION(SIZE(n)) :: arg2d1 - REAL(pr), DIMENSION(SIZE(n)) :: arg2d0 - REAL(pr), DIMENSION(SIZE(n)) :: arg2d0d - REAL(pr), DIMENSION(SIZE(n)) :: arg2d - REAL(pr), DIMENSION(SIZE(n)) :: arg2dd0 - REAL(pr), DIMENSION(SIZE(n)) :: arg2dd - REAL(pr), DIMENSION(SIZE(n)) :: arg2ddd - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10 - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10d - REAL(pr) :: arg11 - REAL(pr) :: arg11d1 - REAL(pr) :: arg11d0 - REAL(pr) :: arg11d0d - REAL(pr) :: arg11d - REAL(pr) :: arg11dd0 - REAL(pr) :: arg11dd - REAL(pr) :: arg11ddd - REAL(pr), DIMENSION(SIZE(n)) :: temp - REAL(pr), DIMENSION(SIZE(n)) :: tempd0 - REAL(pr), DIMENSION(SIZE(n)) :: tempd - REAL(pr), DIMENSION(SIZE(n)) :: tempdd - REAL(pr) :: temp0 - REAL(pr) :: temp0d0 - REAL(pr) :: temp0d - REAL(pr) :: temp0dd - REAL(pr) :: temp1 - REAL(pr) :: temp2 - REAL(pr) :: temp2d0 - REAL(pr) :: temp2d - REAL(pr) :: temp3 - REAL(pr) :: temp3d0 - REAL(pr) :: temp3d - REAL(pr) :: temp3dd - REAL(pr) :: temp4 - REAL(pr) :: temp4d0 - REAL(pr) :: temp4d - REAL(pr) :: temp4dd - REAL(pr) :: temp5 - REAL(pr) :: temp5d0 - REAL(pr) :: temp5d - REAL(pr) :: temp5dd - REAL(pr), DIMENSION(SIZE(n)) :: temp6 - REAL(pr), DIMENSION(SIZE(n)) :: temp6d - REAL(pr) :: temp7 - REAL(pr) :: temp7d - REAL(pr) :: temp8 - REAL(pr) :: temp8d - REAL(pr) :: temp9 - REAL(pr) :: temp9d - REAL(pr) :: temp10 - REAL(pr) :: temp10d - REAL(pr) :: temp11 - REAL(pr) :: temp12 - REAL(pr) :: temp12d - REAL(pr) :: temp13 - REAL(pr) :: temp13d - REAL(pr), DIMENSION(SIZE(n)) :: temp14 - REAL(pr) :: temp15 - REAL(pr) :: temp16 - REAL(pr) :: temp17 - REAL(pr) :: temp18 - REAL(pr) :: temp19 - LOGICAL, DIMENSION(SIZE(n)) :: mask - LOGICAL, DIMENSION(SIZE(n)) :: mask0 - LOGICAL, DIMENSION(SIZE(n)) :: mask1 - LOGICAL, DIMENSION(SIZE(n)) :: mask2 - LOGICAL, DIMENSION(SIZE(n)) :: mask3 - LOGICAL, DIMENSION(SIZE(n)) :: mask4 - LOGICAL, DIMENSION(SIZE(n)) :: mask5 - LOGICAL, DIMENSION(SIZE(n)) :: mask6 - LOGICAL, DIMENSION(SIZE(n)) :: mask7 - LOGICAL, DIMENSION(SIZE(n)) :: mask8 - LOGICAL, DIMENSION(SIZE(n)) :: mask9 - LOGICAL, DIMENSION(SIZE(n)) :: mask10 - LOGICAL, DIMENSION(SIZE(n)) :: mask11 - tc = model%tc - pc = model%pc - w = model%w - ac = model%ac - b = model%b - k = model%k - r = model%r - kij = model%kij - lij = model%lij - arg1d(:) = td/tc - arg1d0(:) = td0/tc - arg1d1(:) = td1/tc - arg1(:) = t/tc - temp14 = SQRT(arg1(:)) - mask(:) = arg1(:) .EQ. 0.0 - WHERE (mask(:)) - temp6d = 0.0_pr - ELSEWHERE - temp6d = arg1d1(:)/(2.0*temp14) - END WHERE - temp6 = temp14 - mask0(:) = arg1(:) .EQ. 0.0 - WHERE (mask0(:)) tempd = 0.0_pr - tempdd = 0.0_pr - mask1(:) = .NOT.arg1(:) .EQ. 0.0 - WHERE (mask1(:)) - temp14 = arg1d0(:)/(2.0*temp6) - tempdd = -(temp14*temp6d/temp6) - tempd = temp14 - END WHERE - tempd0 = temp6d - temp = temp6 - mask2(:) = arg1(:) .EQ. 0.0 - WHERE (mask2(:)) result1d = 0.0_pr - result1dd = 0.0_pr - mask3(:) = .NOT.arg1(:) .EQ. 0.0 - WHERE (mask3(:)) - temp14 = arg1d(:)/(2.0*temp) - temp6d = -(temp14*tempd0/temp) - temp6 = temp14 - END WHERE - result1ddd = 0.0_pr - mask4(:) = .NOT.arg1(:) .EQ. 0.0 - WHERE (mask4(:)) - temp14 = temp6*tempd/temp - result1ddd = -((tempd*temp6d+temp6*tempdd-temp14*tempd0)/temp) - result1dd = -temp14 - END WHERE - result1dd0 = 0.0_pr - mask5(:) = .NOT.arg1(:) .EQ. 0.0 - WHERE (mask5(:)) - result1dd0 = temp6d - result1d = temp6 - END WHERE - result1d0d = tempdd - result1d0 = tempd - result1d1 = tempd0 - result1 = temp - arg2ddd(:) = -(ac*2*k*((k*(1.0-result1)+1.0)*result1ddd-result1dd*k*& -& result1d1-k*(result1d0*result1dd0+result1d*result1d0d))) - arg2dd(:) = -(ac*2*k*((k*(1.0-result1)+1.0)*result1dd-result1d*k*& -& result1d0)) - arg2dd0(:) = -(ac*2*k*((k*(1.0-result1)+1.0)*result1dd0-result1d*k*& -& result1d1)) - arg2d(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d) - arg2d0d(:) = -(ac*2*k*((k*(1.0-result1)+1.0)*result1d0d-result1d0*k*& -& result1d1)) - arg2d0(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d0) - arg2d1(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d1) - arg2(:) = ac*(1.0+k*(1.0-result1))**2 - temp14 = SQRT(arg2(:)) - mask6(:) = arg2(:) .EQ. 0.0 - WHERE (mask6(:)) - temp6d = 0.0_pr - ELSEWHERE - temp6d = arg2d1(:)/(2.0*temp14) - END WHERE - temp6 = temp14 - mask7(:) = arg2(:) .EQ. 0.0 - WHERE (mask7(:)) - tempdd = 0.0_pr - tempd = 0.0_pr - ELSEWHERE - temp14 = arg2d0(:)/(2.0*temp6) - tempdd = (arg2d0d(:)-temp14*2.0*temp6d)/(2.0*temp6) - tempd = temp14 - END WHERE - tempd0 = temp6d - temp = temp6 - mask8(:) = arg2(:) .EQ. 0.0 - WHERE (mask8(:)) ad = 0.0_pr - add = 0.0_pr - mask9(:) = .NOT.arg2(:) .EQ. 0.0 - WHERE (mask9(:)) - temp14 = arg2d(:)/(2.0*temp) - temp6d = (arg2dd0(:)-temp14*2.0*tempd0)/(2.0*temp) - temp6 = temp14 - END WHERE - addd = 0.0_pr - mask10(:) = .NOT.arg2(:) .EQ. 0.0 - WHERE (mask10(:)) - temp14 = (arg2dd(:)-2.0*temp6*tempd)/(2.0*temp) - addd = (arg2ddd(:)-2.0*(tempd*temp6d+temp6*tempdd)-temp14*2.0*& -& tempd0)/(2.0*temp) - add = temp14 - END WHERE - add0 = 0.0_pr - mask11(:) = .NOT.arg2(:) .EQ. 0.0 - WHERE (mask11(:)) - add0 = temp6d - ad = temp6 - END WHERE - ad0d = tempdd - ad0 = tempd - ad1 = tempd0 - a = temp - amix = 0.0 - bmix = 0.0 - bmixd = 0.0_pr - amixd = 0.0_pr - amixd0 = 0.0_pr - amixdd = 0.0_pr - amixddd = 0.0_pr - amixd1 = 0.0_pr - amixd0d = 0.0_pr - amixdd0 = 0.0_pr - DO i=1,SIZE(n)-1 - DO j=i+1,SIZE(n) - nijd = n(j)*nd(i) + n(i)*nd(j) - nij = n(i)*n(j) - temp0 = 2*(-kij(i, j)+1) - temp7d = nijd*ad1(i) + nij*add0(i) - temp7 = nijd*a(i) + nij*ad(i) - temp15 = nijd*ad0(i) + nij*add(i) - amixddd = amixddd + temp0*(ad0(j)*temp7d+temp7*ad0d(j)+temp15*& -& ad1(j)+a(j)*(nijd*ad0d(i)+nij*addd(i))+nij*(ad0(i)*add0(j)+ad(& -& j)*ad0d(i)+add(j)*ad1(i)+a(i)*addd(j))) - amixdd = amixdd + temp0*(temp7*ad0(j)+a(j)*temp15+nij*(ad(j)*ad0& -& (i)+a(i)*add(j))) - amixdd0 = amixdd0 + temp0*(temp7*ad1(j)+a(j)*temp7d+nij*(ad(j)*& -& ad1(i)+a(i)*add0(j))) - amixd = amixd + temp0*(a(j)*temp7+nij*(a(i)*ad(j))) - amixd0d = amixd0d + temp0*nij*(ad0(i)*ad1(j)+a(j)*ad0d(i)+ad0(j)& -& *ad1(i)+a(i)*ad0d(j)) - amixd0 = amixd0 + temp0*nij*(a(j)*ad0(i)+a(i)*ad0(j)) - amixd1 = amixd1 + temp0*nij*(a(j)*ad1(i)+a(i)*ad1(j)) - amix = amix + temp0*(nij*a(i)*a(j)) - bmixd = bmixd + (b(i)+b(j))*(1-lij(i, j))*nijd - bmix = bmix + nij*(b(i)+b(j))*(1-lij(i, j)) - END DO - END DO - arg1ddd(:) = 2**2*n*nd*(ad0*ad1+a*ad0d) + n**2*2*(ad0*add0+ad*ad0d+& -& add*ad1+a*addd) - arg1dd(:) = n*2**2*nd*a*ad0 + n**2*2*(ad*ad0+a*add) - arg1dd0(:) = n*2**2*nd*a*ad1 + n**2*2*(ad*ad1+a*add0) - arg1d(:) = a**2*2*n*nd + n**2*2*a*ad - arg1d0d(:) = n**2*2*(ad0*ad1+a*ad0d) - arg1d0(:) = n**2*2*a*ad0 - arg1d1(:) = n**2*2*a*ad1 - arg1(:) = n**2*a**2 - amixddd = amixddd + SUM(arg1ddd(:)) - amixdd = amixdd + SUM(arg1dd(:)) - amixdd0 = amixdd0 + SUM(arg1dd0(:)) - amixd = amixd + SUM(arg1d(:)) - amixd0d = amixd0d + SUM(arg1d0d(:)) - amixd0 = amixd0 + SUM(arg1d0(:)) - amixd1 = amixd1 + SUM(arg1d1(:)) - amix = amix + SUM(arg1(:)) - arg10d(:) = model%b*2*n*nd - arg10(:) = n**2*model%b - bmixd = bmixd + SUM(arg10d(:)) - bmix = bmix + SUM(arg10(:)) - temp0 = SUM(n) - bmixd = (bmixd-bmix*SUM(nd)/temp0)/temp0 - bmix = bmix/temp0 - temp15 = bmix*vd/v - temp7d = -(temp15*vd1/v) - temp7 = temp15 - temp15 = (bmixd-temp7)/v - temp8d = (-temp7d-temp15*vd1)/v - temp8 = temp15 - temp15 = (temp7/v-temp8)/v - b_vddd = vd0*((temp7d-temp7*vd1/v)/v-temp8d-temp15*vd1)/v - b_vdd = vd0*temp15 - b_vdd0 = temp8d - b_vd = temp8 - temp15 = bmix*vd0/(v*v) - b_vd0d = temp15*2*vd1/v - b_vd0 = -temp15 - b_vd1 = -(bmix*vd1/v**2) - b_v = bmix/v - temp15 = (del1*b_v+1.0)/(del2*b_v+1.0) - temp8d = (del1-temp15*del2)*b_vd1/(del2*b_v+1.0) - temp8 = temp15 - temp15 = b_vd0/(del2*b_v+1.0) - temp0dd = (del1-del2*temp8)*(b_vd0d-temp15*del2*b_vd1)/(del2*b_v+1.0& -& ) - temp15*del2*temp8d - temp0d = (del1-del2*temp8)*temp15 - temp0d0 = temp8d - temp0 = temp8 - temp15 = b_vd/(del2*b_v+1.0) - temp8d = (b_vdd0-temp15*del2*b_vd1)/(del2*b_v+1.0) - temp8 = temp15 - temp15 = b_vdd - del2*temp8*b_vd0 - temp16 = (del1-del2*temp0)/(del2*b_v+1.0) - arg11ddd = temp15*(-(del2*temp0d0)-temp16*del2*b_vd1)/(del2*b_v+1.0)& -& + temp16*(b_vddd-del2*(b_vd0*temp8d+temp8*b_vd0d)) - del2*(temp0d*& -& temp8d+temp8*temp0dd) - arg11dd = temp16*temp15 - del2*(temp8*temp0d) - arg11dd0 = (del1-del2*temp0)*temp8d - temp8*del2*temp0d0 - arg11d = (del1-del2*temp0)*temp8 - arg11d0d = temp0dd - arg11d0 = temp0d - arg11d1 = temp0d0 - arg11 = temp0 - temp16 = b_vd0/(-b_v+1.0) - temp0dd = -((b_vd0d+temp16*b_vd1)/(1.0-b_v)) - temp0d = -temp16 - temp0d0 = -(b_vd1/(1.0-b_v)) - temp0 = LOG(-b_v + 1.0) - temp1 = SUM(n) - temp2d = r*bmix*(del1-del2)*td0 - temp2d0 = r*bmix*(del1-del2)*td1 - temp2 = r*(del1-del2)*t*bmix - temp16 = (amixd0-temp2d*amix/temp2)/temp2 - temp3dd = (amixd0d-temp2d*(amixd1-amix*temp2d0/temp2)/temp2-temp16*& -& temp2d0)/temp2 - temp3d = temp16 - temp3d0 = (amixd1-amix*temp2d0/temp2)/temp2 - temp3 = amix/temp2 - temp4dd = (arg11d0d-arg11d0*arg11d1/arg11)/arg11 - temp4d = arg11d0/arg11 - temp4d0 = arg11d1/arg11 - temp4 = LOG(arg11) - temp5dd = -(temp1*temp0dd) - temp4d*temp3d0 - temp3*temp4dd - temp3d& -& *temp4d0 - temp4*temp3dd - temp5d = -(temp1*temp0d) - temp3*temp4d - temp4*temp3d - temp5d0 = -(temp1*temp0d0) - temp3*temp4d0 - temp4*temp3d0 - temp5 = -(temp1*temp0) - temp4*temp3 - temp8d = (temp4d0-temp4*temp2d0/temp2)/temp2 - temp8 = temp4/temp2 - temp7d = bmixd*td1 - temp7 = bmix*td + bmixd*t - temp9d = amixdd0 - r*(del1-del2)*(temp7*temp3d0+temp3*temp7d) - temp9 = amixd - r*(del1-del2)*temp3*temp7 - temp16 = temp3*arg11d/arg11 - temp10d = (arg11d*temp3d0+temp3*arg11dd0-temp16*arg11d1)/arg11 - temp10 = temp16 - temp11 = SUM(nd) - temp16 = b_vd/(-b_v+1.0) - temp12d = (b_vdd0+temp16*b_vd1)/(1.0-b_v) - temp12 = temp16 - temp13d = temp1*temp12d - temp11*temp0d0 - temp10d - temp8*temp9d - & -& temp9*temp8d - temp13 = temp1*temp12 - temp11*temp0 - temp10 - temp9*temp8 - temp16 = temp9/temp2 - temp15 = amixdd - r*(del1-del2)*(temp7*temp3d+bmixd*td0*temp3) - temp17 = (arg11d*temp3d+temp3*arg11dd-temp10*arg11d0)/arg11 - temp18 = (b_vdd+temp12*b_vd0)/(-b_v+1.0) - temp19 = temp1*temp18 - temp11*temp0d - temp17 - temp8*temp15 - (& -& temp4d-temp2d*temp8)*temp16 - arvalddd = model%r*(td0*temp13d+temp19*td1+t*(temp1*(b_vddd+b_vd0*& -& temp12d+temp12*b_vd0d+temp18*b_vd1)/(1.0-b_v)-temp11*temp0dd-(& -& temp3d*arg11dd0+arg11d*temp3dd+arg11dd*temp3d0+temp3*arg11ddd-& -& arg11d0*temp10d-temp10*arg11d0d-temp17*arg11d1)/arg11-temp15*& -& temp8d-temp8*(amixddd-r*(del1-del2)*(temp3d*temp7d+temp7*temp3dd+& -& bmixd*td0*temp3d0))-temp16*(temp4dd-temp2d*temp8d)-(temp4d-temp2d*& -& temp8)*(temp9d-temp16*temp2d0)/temp2)+td*temp5dd) - arvaldd = model%r*(td0*temp13+t*temp19+td*temp5d) - arvaldd0 = model%r*(temp13*td1+t*temp13d+td*temp5d0) - arvald = model%r*(t*temp13+td*temp5) - arvald0d = model%r*(temp5d*td1+t*temp5dd+td0*temp5d0) - arvald0 = model%r*(t*temp5d+temp5*td0) - arvald1 = model%r*(t*temp5d0+temp5*td1) - arval = model%r*(temp5*t) - END SUBROUTINE AR_D_D_D - -! Differentiation of ar_d in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval arvald -! with respect to varying inputs: t v -! RW status of diff variables: t:in v:in arval:out arvald:out -! Differentiation of ar in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval -! with respect to varying inputs: n t v -! RW status of diff variables: n:in t:in v:in arval:out - SUBROUTINE AR_D_D(model, n, nd, v, vd0, vd, t, td0, td, arval, arvald0& -& , arvald, arvaldd) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:), v, t - REAL(pr), INTENT(IN) :: vd0, td0 - REAL(pr), INTENT(IN) :: nd(:), vd, td - REAL(pr), INTENT(OUT) :: arval - REAL(pr), INTENT(OUT) :: arvald0 - REAL(pr), INTENT(OUT) :: arvald - REAL(pr), INTENT(OUT) :: arvaldd - REAL(pr) :: amix, a(SIZE(n)), ai(SIZE(n)), z2(SIZE(n)), nij - REAL(pr) :: amixd0, ad0(SIZE(n)) - REAL(pr) :: amixd, ad(SIZE(n)), nijd - REAL(pr) :: amixdd, add(SIZE(n)) - REAL(pr) :: bmix - REAL(pr) :: bmixd - REAL(pr) :: b_v - REAL(pr) :: b_vd0 - REAL(pr) :: b_vd - REAL(pr) :: b_vdd - REAL(pr) :: aij(SIZE(n), SIZE(n)), bij(SIZE(n), SIZE(n)) - INTEGER :: i, j - REAL(pr) :: tc(SIZE(n)), pc(SIZE(n)), w(SIZE(n)) - REAL(pr) :: ac(SIZE(n)), b(SIZE(n)), k(SIZE(n)), r - REAL(pr) :: del1, del2 - REAL(pr) :: kij(SIZE(n), SIZE(n)), lij(SIZE(n), SIZE(n)) - INTRINSIC SQRT - INTRINSIC SUM - INTRINSIC LOG - INTRINSIC SIZE - REAL(pr), DIMENSION(SIZE(n)) :: arg1 - REAL(pr), DIMENSION(SIZE(n)) :: arg1d0 - REAL(pr), DIMENSION(SIZE(n)) :: arg1d - REAL(pr), DIMENSION(SIZE(n)) :: arg1dd - REAL(pr), DIMENSION(SIZE(n)) :: result1 - REAL(pr), DIMENSION(SIZE(n)) :: result1d0 - REAL(pr), DIMENSION(SIZE(n)) :: result1d - REAL(pr), DIMENSION(SIZE(n)) :: result1dd - REAL(pr), DIMENSION(SIZE(n)) :: arg2 - REAL(pr), DIMENSION(SIZE(n)) :: arg2d0 - REAL(pr), DIMENSION(SIZE(n)) :: arg2d - REAL(pr), DIMENSION(SIZE(n)) :: arg2dd - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10 - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10d - REAL(pr) :: arg11 - REAL(pr) :: arg11d0 - REAL(pr) :: arg11d - REAL(pr) :: arg11dd - REAL(pr), DIMENSION(SIZE(n)) :: temp - REAL(pr), DIMENSION(SIZE(n)) :: tempd - REAL(pr) :: temp0 - REAL(pr) :: temp0d - REAL(pr) :: temp1 - REAL(pr) :: temp2 - REAL(pr) :: temp2d - REAL(pr) :: temp3 - REAL(pr) :: temp3d - REAL(pr) :: temp4 - REAL(pr) :: temp4d - REAL(pr) :: temp5 - REAL(pr) :: temp5d - REAL(pr), DIMENSION(SIZE(n)) :: temp6 - REAL(pr) :: temp7 - REAL(pr) :: temp8 - REAL(pr) :: temp9 - REAL(pr) :: temp10 - REAL(pr) :: temp11 - REAL(pr) :: temp12 - REAL(pr) :: temp13 - LOGICAL, DIMENSION(SIZE(n)) :: mask - LOGICAL, DIMENSION(SIZE(n)) :: mask0 - LOGICAL, DIMENSION(SIZE(n)) :: mask1 - LOGICAL, DIMENSION(SIZE(n)) :: mask2 - LOGICAL, DIMENSION(SIZE(n)) :: mask3 - LOGICAL, DIMENSION(SIZE(n)) :: mask4 - tc = model%tc - pc = model%pc - w = model%w - ac = model%ac - b = model%b - k = model%k - r = model%r - kij = model%kij - lij = model%lij - arg1d(:) = td/tc - arg1d0(:) = td0/tc - arg1(:) = t/tc - temp6 = SQRT(arg1(:)) - mask(:) = arg1(:) .EQ. 0.0 - WHERE (mask(:)) - tempd = 0.0_pr - ELSEWHERE - tempd = arg1d0(:)/(2.0*temp6) - END WHERE - temp = temp6 - mask0(:) = arg1(:) .EQ. 0.0 - WHERE (mask0(:)) result1d = 0.0_pr - result1dd = 0.0_pr - mask1(:) = .NOT.arg1(:) .EQ. 0.0 - WHERE (mask1(:)) - temp6 = arg1d(:)/(2.0*temp) - result1dd = -(temp6*tempd/temp) - result1d = temp6 - END WHERE - result1d0 = tempd - result1 = temp - arg2dd(:) = -(ac*2*k*((k*(1.0-result1)+1.0)*result1dd-result1d*k*& -& result1d0)) - arg2d(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d) - arg2d0(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d0) - arg2(:) = ac*(1.0+k*(1.0-result1))**2 - temp6 = SQRT(arg2(:)) - mask2(:) = arg2(:) .EQ. 0.0 - WHERE (mask2(:)) - tempd = 0.0_pr - ELSEWHERE - tempd = arg2d0(:)/(2.0*temp6) - END WHERE - temp = temp6 - mask3(:) = arg2(:) .EQ. 0.0 - WHERE (mask3(:)) ad = 0.0_pr - add = 0.0_pr - mask4(:) = .NOT.arg2(:) .EQ. 0.0 - WHERE (mask4(:)) - temp6 = arg2d(:)/(2.0*temp) - add = (arg2dd(:)-temp6*2.0*tempd)/(2.0*temp) - ad = temp6 - END WHERE - ad0 = tempd - a = temp - amix = 0.0 - bmix = 0.0 - bmixd = 0.0_pr - amixd = 0.0_pr - amixd0 = 0.0_pr - amixdd = 0.0_pr - DO i=1,SIZE(n)-1 - DO j=i+1,SIZE(n) - nijd = n(j)*nd(i) + n(i)*nd(j) - nij = n(i)*n(j) - temp0 = 2*(-kij(i, j)+1) - temp7 = nijd*a(i) + nij*ad(i) - amixdd = amixdd + temp0*(temp7*ad0(j)+a(j)*(nijd*ad0(i)+nij*add(& -& i))+nij*(ad(j)*ad0(i)+a(i)*add(j))) - amixd = amixd + temp0*(a(j)*temp7+nij*(a(i)*ad(j))) - amixd0 = amixd0 + temp0*nij*(a(j)*ad0(i)+a(i)*ad0(j)) - amix = amix + temp0*(nij*a(i)*a(j)) - bmixd = bmixd + (b(i)+b(j))*(1-lij(i, j))*nijd - bmix = bmix + nij*(b(i)+b(j))*(1-lij(i, j)) - END DO - END DO - arg1dd(:) = n*2**2*nd*a*ad0 + n**2*2*(ad*ad0+a*add) - arg1d(:) = a**2*2*n*nd + n**2*2*a*ad - arg1d0(:) = n**2*2*a*ad0 - arg1(:) = n**2*a**2 - amixdd = amixdd + SUM(arg1dd(:)) - amixd = amixd + SUM(arg1d(:)) - amixd0 = amixd0 + SUM(arg1d0(:)) - amix = amix + SUM(arg1(:)) - arg10d(:) = model%b*2*n*nd - arg10(:) = n**2*model%b - bmixd = bmixd + SUM(arg10d(:)) - bmix = bmix + SUM(arg10(:)) - temp0 = SUM(n) - bmixd = (bmixd-bmix*SUM(nd)/temp0)/temp0 - bmix = bmix/temp0 - temp7 = bmix*vd/v - temp8 = (bmixd-temp7)/v - b_vdd = (temp7/v-temp8)*vd0/v - b_vd = temp8 - b_vd0 = -(bmix*vd0/v**2) - b_v = bmix/v - temp8 = (del1*b_v+1.0)/(del2*b_v+1.0) - temp0d = (del1-temp8*del2)*b_vd0/(del2*b_v+1.0) - temp0 = temp8 - temp8 = b_vd/(del2*b_v+1.0) - arg11dd = (del1-del2*temp0)*(b_vdd-temp8*del2*b_vd0)/(del2*b_v+1.0) & -& - temp8*del2*temp0d - arg11d = (del1-del2*temp0)*temp8 - arg11d0 = temp0d - arg11 = temp0 - temp0d = -(b_vd0/(1.0-b_v)) - temp0 = LOG(-b_v + 1.0) - temp1 = SUM(n) - temp2d = r*bmix*(del1-del2)*td0 - temp2 = r*(del1-del2)*t*bmix - temp3d = (amixd0-amix*temp2d/temp2)/temp2 - temp3 = amix/temp2 - temp4d = arg11d0/arg11 - temp4 = LOG(arg11) - temp5d = -(temp1*temp0d) - temp3*temp4d - temp4*temp3d - temp5 = -(temp1*temp0) - temp4*temp3 - temp8 = temp4/temp2 - temp7 = bmix*td + bmixd*t - temp9 = amixd - r*(del1-del2)*temp3*temp7 - temp10 = temp3*arg11d/arg11 - temp11 = SUM(nd) - temp12 = b_vd/(-b_v+1.0) - temp13 = temp1*temp12 - temp11*temp0 - temp10 - temp9*temp8 - arvaldd = model%r*(temp13*td0+t*(temp1*(b_vdd+temp12*b_vd0)/(1.0-b_v& -& )-temp11*temp0d-(arg11d*temp3d+temp3*arg11dd-temp10*arg11d0)/arg11& -& -temp8*(amixdd-r*(del1-del2)*(temp7*temp3d+temp3*bmixd*td0))-temp9& -& *(temp4d-temp8*temp2d)/temp2)+td*temp5d) - arvald = model%r*(t*temp13+td*temp5) - arvald0 = model%r*(t*temp5d+temp5*td0) - arval = model%r*(temp5*t) - END SUBROUTINE AR_D_D - -! Differentiation of ar_d in reverse (adjoint) mode (with options noISIZE): -! gradient of useful results: nd n t v arval arvald vd td -! with respect to varying inputs: nd n t v arval arvald vd td -! RW status of diff variables: nd:incr n:incr t:incr v:incr arval:in-zero -! arvald:in-zero vd:incr td:incr -! Differentiation of ar in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval -! with respect to varying inputs: n t v -! RW status of diff variables: n:in t:in v:in arval:out - SUBROUTINE AR_D_B(model, n, nb, nd, ndb, v, vb, vd, vdb, t, tb, td, & -& tdb, arval, arvalb, arvald, arvaldb) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:), v, t - REAL(pr) :: nb(:), vb, tb - REAL(pr), INTENT(IN) :: nd(:), vd, td - REAL(pr) :: ndb(:), vdb, tdb - REAL(pr) :: arval - REAL(pr) :: arvalb - REAL(pr) :: arvald - REAL(pr) :: arvaldb - REAL(pr) :: amix, a(SIZE(n)), ai(SIZE(n)), z2(SIZE(n)), nij - REAL(pr) :: amixb, ab(SIZE(n)), nijb - REAL(pr) :: amixd, ad(SIZE(n)), nijd - REAL(pr) :: amixdb, adb(SIZE(n)), nijdb - REAL(pr) :: bmix - REAL(pr) :: bmixb - REAL(pr) :: bmixd - REAL(pr) :: bmixdb - REAL(pr) :: b_v - REAL(pr) :: b_vb - REAL(pr) :: b_vd - REAL(pr) :: b_vdb - REAL(pr) :: aij(SIZE(n), SIZE(n)), bij(SIZE(n), SIZE(n)) - INTEGER :: i, j - REAL(pr) :: tc(SIZE(n)), pc(SIZE(n)), w(SIZE(n)) - REAL(pr) :: ac(SIZE(n)), b(SIZE(n)), k(SIZE(n)), r - REAL(pr) :: del1, del2 - REAL(pr) :: kij(SIZE(n), SIZE(n)), lij(SIZE(n), SIZE(n)) - INTRINSIC SQRT - INTRINSIC SUM - INTRINSIC LOG - INTRINSIC SIZE - REAL(pr), DIMENSION(SIZE(n)) :: arg1 - REAL(pr), DIMENSION(SIZE(n)) :: arg1b - REAL(pr), DIMENSION(SIZE(n)) :: arg1d - REAL(pr), DIMENSION(SIZE(n)) :: arg1db - REAL(pr), DIMENSION(SIZE(n)) :: result1 - REAL(pr), DIMENSION(SIZE(n)) :: result1b - REAL(pr), DIMENSION(SIZE(n)) :: result1d - REAL(pr), DIMENSION(SIZE(n)) :: result1db - REAL(pr), DIMENSION(SIZE(n)) :: arg2 - REAL(pr), DIMENSION(SIZE(n)) :: arg2b - REAL(pr), DIMENSION(SIZE(n)) :: arg2d - REAL(pr), DIMENSION(SIZE(n)) :: arg2db - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10 - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10b - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10d - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10db - REAL(pr) :: arg11 - REAL(pr) :: arg11b - REAL(pr) :: arg11d - REAL(pr) :: arg11db - REAL(pr), DIMENSION(SIZE(n)) :: temp - REAL(pr), DIMENSION(SIZE(n)) :: tempb - REAL(pr) :: temp0 - REAL(pr) :: temp0b - REAL(pr) :: temp1 - REAL(pr) :: temp1b - REAL(pr) :: temp2 - REAL(pr) :: temp2b - REAL(pr) :: temp3 - REAL(pr) :: temp3b - REAL(pr) :: temp4 - REAL(pr) :: temp4b - REAL(pr) :: temp5 - REAL(pr) :: temp5b - LOGICAL, DIMENSION(SIZE(n)) :: mask - LOGICAL, DIMENSION(SIZE(n)) :: mask0 - REAL(pr), DIMENSION(SIZE(n)) :: tempb0 - REAL(pr) :: temp6 - REAL(pr) :: tempb1 - REAL(pr) :: temp7 - REAL(pr) :: tempb2 - REAL(pr), DIMENSION(SIZE(n, 1)) :: tempb3 - REAL(pr) :: temp8 - REAL(pr) :: tempb4 - REAL(pr) :: temp9 - REAL(pr) :: temp10 - REAL(pr) :: temp11 - REAL(pr) :: temp12 - REAL(pr) :: tempb5 - REAL(pr) :: tempb6 - REAL(pr) :: tempb7 - REAL(pr) :: tempb8 - INTEGER :: ad_from - INTEGER :: ad_to - INTEGER :: ad_to0 - INTEGER :: arg12 - LOGICAL, DIMENSION(SIZE(n)) :: mask1 - LOGICAL, DIMENSION(SIZE(n)) :: mask2 - LOGICAL, DIMENSION(SIZE(n)) :: mask3 - LOGICAL, DIMENSION(SIZE(n)) :: mask4 - tc = model%tc - ac = model%ac - b = model%b - k = model%k - r = model%r - kij = model%kij - lij = model%lij - arg1d(:) = td/tc - arg1(:) = t/tc - temp = SQRT(arg1(:)) - mask(:) = arg1(:) .EQ. 0.0 - WHERE (mask(:)) - result1d = 0.0_pr - ELSEWHERE - result1d = arg1d(:)/(2.0*temp) - END WHERE - result1 = temp - arg2d(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d) - arg2(:) = ac*(1.0+k*(1.0-result1))**2 - arg12 = SIZE(n) - CALL PUSHREAL8ARRAY(temp, arg12) - temp = SQRT(arg2(:)) - mask0(:) = arg2(:) .EQ. 0.0 - WHERE (mask0(:)) - ad = 0.0_pr - ELSEWHERE - ad = arg2d(:)/(2.0*temp) - END WHERE - a = temp - amix = 0.0 - bmix = 0.0 - bmixd = 0.0_pr - amixd = 0.0_pr - DO i=1,SIZE(n)-1 - ad_from = i + 1 - DO j=ad_from,SIZE(n) - nijd = n(j)*nd(i) + n(i)*nd(j) - nij = n(i)*n(j) - temp0 = 2*(-kij(i, j)+1) - amixd = amixd + temp0*(a(j)*(a(i)*nijd+nij*ad(i))+nij*a(i)*ad(j)& -& ) - amix = amix + temp0*(nij*a(i)*a(j)) - bmixd = bmixd + (b(i)+b(j))*(1-lij(i, j))*nijd - bmix = bmix + nij*(b(i)+b(j))*(1-lij(i, j)) - END DO - CALL PUSHINTEGER4(j - 1) - CALL PUSHINTEGER4(ad_from) - END DO - CALL PUSHINTEGER4(i - 1) - arg1d(:) = a**2*2*n*nd + n**2*2*a*ad - arg1(:) = n**2*a**2 - amixd = amixd + SUM(arg1d(:)) - amix = amix + SUM(arg1(:)) - arg10d(:) = model%b*2*n*nd - arg10(:) = n**2*model%b - bmixd = bmixd + SUM(arg10d(:)) - bmix = bmix + SUM(arg10(:)) - temp0 = SUM(n) - CALL PUSHREAL8(bmixd) - bmixd = (bmixd-bmix*SUM(nd)/temp0)/temp0 - CALL PUSHREAL8(bmix) - bmix = bmix/temp0 - b_vd = (bmixd-bmix*vd/v)/v - b_v = bmix/v - CALL PUSHREAL8(temp0) - temp0 = (del1*b_v+1.0)/(del2*b_v+1.0) - arg11d = (del1-temp0*del2)*b_vd/(del2*b_v+1.0) - arg11 = temp0 - temp0 = LOG(-b_v + 1.0) - temp1 = SUM(n) - temp2 = r*(del1-del2)*t*bmix - temp3 = amix/temp2 - temp4 = LOG(arg11) - temp5 = -(temp1*temp0) - temp4*temp3 - tempb4 = model%r*arvaldb - temp12 = temp1*b_vd/(-b_v+1.0) - temp11 = SUM(nd) - temp10 = temp3*arg11d/arg11 - temp6 = bmix*td + t*bmixd - temp9 = amixd - r*(del1-del2)*temp3*temp6 - temp7 = temp4/temp2 - temp5b = t*model%r*arvalb + td*tempb4 - tb = tb + temp5*model%r*arvalb + (temp12-temp0*temp11-temp10-temp9*& -& temp7)*tempb4 - tempb5 = t*tempb4 - tempb6 = tempb5/(1.0-b_v) - temp0b = -(temp11*tempb5) - temp1*temp5b - ndb = ndb - temp0*tempb5 - tempb7 = -(tempb5/arg11) - amixdb = -(temp7*tempb5) - tempb8 = r*(del1-del2)*temp7*tempb5 - tempb2 = -(temp9*tempb5/temp2) - temp4b = tempb2 - temp3*temp5b - temp3b = temp6*tempb8 + arg11d*tempb7 - temp4*temp5b - temp2b = -(temp7*tempb2) - amix*temp3b/temp2**2 - tempb1 = temp3*tempb8 - tdb = tdb + temp5*tempb4 + bmix*tempb1 - arg11db = temp3*tempb7 - arg11b = temp4b/arg11 - temp10*tempb7 - temp1b = b_vd*tempb6 - temp0*temp5b - amixb = temp3b/temp2 - tempb4 = r*(del1-del2)*temp2b - bmixb = td*tempb1 + t*tempb4 - tb = tb + bmixd*tempb1 + bmix*tempb4 - temp0 = (del1*b_v+1.0)/(del2*b_v+1.0) - temp8 = b_vd/(del2*b_v+1.0) - tempb4 = (del1-del2*temp0)*arg11db/(del2*b_v+1.0) - b_vdb = temp1*tempb6 + tempb4 - b_vb = temp12*tempb6 - temp0b/(1.0-b_v) - del2*temp8*tempb4 - temp0b = arg11b - del2*temp8*arg11db - CALL POPREAL8(temp0) - tempb4 = temp0b/(del2*b_v+1.0) - b_vb = b_vb + (del1-del2*(del1*b_v+1.0)/(del2*b_v+1.0))*tempb4 - temp7 = bmix*vd/v - tempb4 = b_vdb/v - bmixdb = t*tempb1 + tempb4 - tempb2 = -(tempb4/v) - bmixb = bmixb + b_vb/v + vd*tempb2 - vb = vb - bmix*b_vb/v**2 - (bmixd-temp7)*tempb4/v - temp7*tempb2 - vdb = vdb + bmix*tempb2 - CALL POPREAL8(bmix) - CALL POPREAL8(bmixd) - temp6 = bmix/temp0 - temp8 = SUM(nd) - tempb2 = bmixdb/temp0 - bmixdb = tempb2 - tempb1 = -(temp8*tempb2/temp0) - temp0b = -(bmix*bmixb/temp0**2) - (bmixd-temp8*temp6)*tempb2/temp0 -& -& temp6*tempb1 - bmixb = bmixb/temp0 + tempb1 - arg10b = 0.0_pr - arg10b = bmixb - arg10db = 0.0_pr - arg10db = bmixdb - arg1b = 0.0_pr - arg1b = amixb - arg1db = 0.0_pr - arg1db = amixdb - ab = 0.0_pr - arg1(:) = t/tc - adb = 0.0_pr - tempb3 = 2*a**2*arg1db - nb = nb + temp1b + temp0b + 2*n*model%b*arg10b + nd*model%b*2*& -& arg10db + 2*n*a**2*arg1b + 2**2*n*a*ad*arg1db + nd*tempb3 - ndb = ndb + n*model%b*2*arg10db - temp6*tempb2 + n*tempb3 - tempb0 = 2*n**2*arg1db - ab = 2*a*n**2*arg1b + 2**2*a*n*nd*arg1db + ad*tempb0 - adb = a*tempb0 - CALL POPINTEGER4(ad_to0) - DO i=ad_to0,1,-1 - CALL POPINTEGER4(ad_from) - CALL POPINTEGER4(ad_to) - DO j=ad_to,ad_from,-1 - temp0 = 2*(-kij(i, j)+1) - tempb2 = a(j)*temp0*amixb - nij = n(i)*n(j) - nijb = (1-lij(i, j))*(b(i)+b(j))*bmixb + a(i)*tempb2 - nijd = n(j)*nd(i) + n(i)*nd(j) - ab(j) = ab(j) + nij*a(i)*temp0*amixb - ab(i) = ab(i) + nij*tempb2 - tempb1 = temp0*amixdb - ab(j) = ab(j) + (a(i)*nijd+nij*ad(i))*tempb1 - tempb2 = a(j)*tempb1 - nijdb = (1-lij(i, j))*(b(i)+b(j))*bmixdb + a(i)*tempb2 - nijb = nijb + a(i)*ad(j)*tempb1 + ad(i)*tempb2 - ab(i) = ab(i) + nij*ad(j)*tempb1 + nijd*tempb2 - adb(j) = adb(j) + nij*a(i)*tempb1 - adb(i) = adb(i) + nij*tempb2 - nb(i) = nb(i) + n(j)*nijb - nb(j) = nb(j) + n(i)*nijb + nd(i)*nijdb - ndb(i) = ndb(i) + n(j)*nijdb - nb(i) = nb(i) + nd(j)*nijdb - ndb(j) = ndb(j) + n(i)*nijdb - END DO - END DO - tempb = 0.0_pr - tempb = ab - arg2db = 0.0_pr - arg2b = 0.0_pr - result1b = 0.0_pr - result1db = 0.0_pr - arg1d(:) = td/tc - arg1db = 0.0_pr - arg1b = 0.0_pr - mask1(:) = .NOT.mask0(:) - WHERE (mask1(:)) - tempb0 = adb/(2.0*temp) - arg2db = tempb0 - tempb = tempb - arg2d*tempb0/temp - END WHERE - tempb0 = -(ac*2*k*arg2db) - arg12 = SIZE(n) - CALL POPREAL8ARRAY(temp, arg12) - mask2 = arg2 .EQ. 0.0 - WHERE (mask2) - arg2b = 0.0_pr - ELSEWHERE - arg2b = tempb/(2.0*SQRT(arg2)) - END WHERE - result1b = -(k*2*(k*(1.0-result1)+1.0)*ac*arg2b) - k*result1d*tempb0 - result1db = (k*(1.0-result1)+1.0)*tempb0 - tempb = 0.0_pr - tempb = result1b - mask3(:) = .NOT.mask(:) - WHERE (mask3(:)) - tempb0 = result1db/(2.0*temp) - arg1db = tempb0 - tempb = tempb - arg1d*tempb0/temp - END WHERE - mask4 = arg1 .EQ. 0.0 - WHERE (mask4) - arg1b = 0.0_pr - ELSEWHERE - arg1b = tempb/(2.0*SQRT(arg1)) - END WHERE - tb = tb + SUM(arg1b/tc) - tdb = tdb + SUM(arg1db/tc) - arvalb = 0.0_pr - arvaldb = 0.0_pr - END SUBROUTINE AR_D_B - -! Differentiation of ar in forward (tangent) mode (with options noISIZE): -! variations of useful results: arval -! with respect to varying inputs: n t v -! RW status of diff variables: n:in t:in v:in arval:out - SUBROUTINE AR_D(model, n, nd, v, vd, t, td, arval, arvald) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:), v, t - REAL(pr), INTENT(IN) :: nd(:), vd, td - REAL(pr), INTENT(OUT) :: arval - REAL(pr), INTENT(OUT) :: arvald - REAL(pr) :: amix, a(SIZE(n)), ai(SIZE(n)), z2(SIZE(n)), nij - REAL(pr) :: amixd, ad(SIZE(n)), nijd - REAL(pr) :: bmix - REAL(pr) :: bmixd - REAL(pr) :: b_v - REAL(pr) :: b_vd - REAL(pr) :: aij(SIZE(n), SIZE(n)), bij(SIZE(n), SIZE(n)) - INTEGER :: i, j - REAL(pr) :: tc(SIZE(n)), pc(SIZE(n)), w(SIZE(n)) - REAL(pr) :: ac(SIZE(n)), b(SIZE(n)), k(SIZE(n)), r - REAL(pr) :: del1, del2 - REAL(pr) :: kij(SIZE(n), SIZE(n)), lij(SIZE(n), SIZE(n)) - INTRINSIC SQRT - INTRINSIC SUM - INTRINSIC LOG - INTRINSIC SIZE - REAL(pr), DIMENSION(SIZE(n)) :: arg1 - REAL(pr), DIMENSION(SIZE(n)) :: arg1d - REAL(pr), DIMENSION(SIZE(n)) :: result1 - REAL(pr), DIMENSION(SIZE(n)) :: result1d - REAL(pr), DIMENSION(SIZE(n)) :: arg2 - REAL(pr), DIMENSION(SIZE(n)) :: arg2d - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10 - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10d - REAL(pr) :: arg11 - REAL(pr) :: arg11d - REAL(pr), DIMENSION(SIZE(n)) :: temp - REAL(pr) :: temp0 - REAL(pr) :: temp1 - REAL(pr) :: temp2 - REAL(pr) :: temp3 - REAL(pr) :: temp4 - REAL(pr) :: temp5 - LOGICAL, DIMENSION(SIZE(n)) :: mask - LOGICAL, DIMENSION(SIZE(n)) :: mask0 - tc = model%tc - pc = model%pc - w = model%w - ac = model%ac - b = model%b - k = model%k - r = model%r - kij = model%kij - lij = model%lij - arg1d(:) = td/tc - arg1(:) = t/tc - temp = SQRT(arg1(:)) - mask(:) = arg1(:) .EQ. 0.0 - WHERE (mask(:)) - result1d = 0.0_pr - ELSEWHERE - result1d = arg1d(:)/(2.0*temp) - END WHERE - result1 = temp - arg2d(:) = -(ac*2*(k*(1.0-result1)+1.0)*k*result1d) - arg2(:) = ac*(1.0+k*(1.0-result1))**2 - temp = SQRT(arg2(:)) - mask0(:) = arg2(:) .EQ. 0.0 - WHERE (mask0(:)) - ad = 0.0_pr - ELSEWHERE - ad = arg2d(:)/(2.0*temp) - END WHERE - a = temp - amix = 0.0 - bmix = 0.0 - bmixd = 0.0_pr - amixd = 0.0_pr - DO i=1,SIZE(n)-1 - DO j=i+1,SIZE(n) - nijd = n(j)*nd(i) + n(i)*nd(j) - nij = n(i)*n(j) - temp0 = 2*(-kij(i, j)+1) - amixd = amixd + temp0*(a(j)*(a(i)*nijd+nij*ad(i))+nij*a(i)*ad(j)& -& ) - amix = amix + temp0*(nij*a(i)*a(j)) - bmixd = bmixd + (b(i)+b(j))*(1-lij(i, j))*nijd - bmix = bmix + nij*(b(i)+b(j))*(1-lij(i, j)) - END DO - END DO - arg1d(:) = a**2*2*n*nd + n**2*2*a*ad - arg1(:) = n**2*a**2 - amixd = amixd + SUM(arg1d(:)) - amix = amix + SUM(arg1(:)) - arg10d(:) = model%b*2*n*nd - arg10(:) = n**2*model%b - bmixd = bmixd + SUM(arg10d(:)) - bmix = bmix + SUM(arg10(:)) - temp0 = SUM(n) - bmixd = (bmixd-bmix*SUM(nd)/temp0)/temp0 - bmix = bmix/temp0 - b_vd = (bmixd-bmix*vd/v)/v - b_v = bmix/v - temp0 = (del1*b_v+1.0)/(del2*b_v+1.0) - arg11d = (del1-temp0*del2)*b_vd/(del2*b_v+1.0) - arg11 = temp0 - temp0 = LOG(-b_v + 1.0) - temp1 = SUM(n) - temp2 = r*(del1-del2)*t*bmix - temp3 = amix/temp2 - temp4 = LOG(arg11) - temp5 = -(temp1*temp0) - temp4*temp3 - arvald = model%r*(t*(temp1*b_vd/(1.0-b_v)-temp0*SUM(nd)-temp3*arg11d& -& /arg11-temp4*(amixd-temp3*r*(del1-del2)*(bmix*td+t*bmixd))/temp2)+& -& temp5*td) - arval = model%r*(temp5*t) - END SUBROUTINE AR_D - -! Differentiation of ar in reverse (adjoint) mode (with options noISIZE): -! gradient of useful results: arval -! with respect to varying inputs: n t v arval -! RW status of diff variables: n:out t:out v:out arval:in-zero - SUBROUTINE AR_B(model, n, nb, v, vb, t, tb, arval, arvalb) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:), v, t - REAL(pr) :: nb(:), vb, tb - REAL(pr) :: arval - REAL(pr) :: arvalb - REAL(pr) :: amix, a(SIZE(n)), ai(SIZE(n)), z2(SIZE(n)), nij - REAL(pr) :: amixb, ab(SIZE(n)), nijb - REAL(pr) :: bmix - REAL(pr) :: bmixb - REAL(pr) :: b_v - REAL(pr) :: b_vb - REAL(pr) :: aij(SIZE(n), SIZE(n)), bij(SIZE(n), SIZE(n)) - INTEGER :: i, j - REAL(pr) :: tc(SIZE(n)), pc(SIZE(n)), w(SIZE(n)) - REAL(pr) :: ac(SIZE(n)), b(SIZE(n)), k(SIZE(n)), r - REAL(pr) :: del1, del2 - REAL(pr) :: kij(SIZE(n), SIZE(n)), lij(SIZE(n), SIZE(n)) - INTRINSIC SQRT - INTRINSIC SUM - INTRINSIC LOG - INTRINSIC SIZE - REAL(pr), DIMENSION(SIZE(n)) :: arg1 - REAL(pr), DIMENSION(SIZE(n)) :: arg1b - REAL(pr), DIMENSION(SIZE(n)) :: result1 - REAL(pr), DIMENSION(SIZE(n)) :: result1b - REAL(pr), DIMENSION(SIZE(n)) :: arg2 - REAL(pr), DIMENSION(SIZE(n)) :: arg2b - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10 - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10b - REAL(pr) :: arg11 - REAL(pr) :: arg11b - REAL(pr) :: temp - REAL(pr) :: tempb - REAL(pr) :: temp0 - REAL(pr) :: temp1 - REAL(pr) :: temp2 - REAL(pr) :: temp3 - REAL(pr) :: temp4 - REAL(pr) :: tempb0 - REAL(pr) :: tempb1 - INTEGER :: ad_from - INTEGER :: ad_to - INTEGER :: ad_to0 - tc = model%tc - ac = model%ac - b = model%b - k = model%k - r = model%r - kij = model%kij - lij = model%lij - arg1(:) = t/tc - result1 = SQRT(arg1(:)) - arg2(:) = ac*(1.0+k*(1.0-result1))**2 - a = SQRT(arg2(:)) - amix = 0.0 - bmix = 0.0 - DO i=1,SIZE(n)-1 - ad_from = i + 1 - DO j=ad_from,SIZE(n) - nij = n(i)*n(j) - amix = amix + 2*nij*(a(i)*a(j))*(1-kij(i, j)) - bmix = bmix + nij*(b(i)+b(j))*(1-lij(i, j)) - END DO - CALL PUSHINTEGER4(j - 1) - CALL PUSHINTEGER4(ad_from) - END DO - CALL PUSHINTEGER4(i - 1) - arg1(:) = n**2*a**2 - amix = amix + SUM(arg1(:)) - arg10(:) = n**2*model%b - bmix = bmix + SUM(arg10(:)) - CALL PUSHREAL8(bmix) - bmix = bmix/SUM(n) - b_v = bmix/v - arg11 = (1.0+del1*b_v)/(1.0+del2*b_v) - nb = 0.0_pr - temp0 = LOG(-b_v + 1.0) - temp1 = SUM(n) - temp2 = r*(del1-del2)*t*bmix - temp3 = amix/temp2 - temp4 = LOG(arg11) - tempb = t*model%r*arvalb - nb = -(temp0*tempb) - b_vb = temp1*tempb/(1.0-b_v) - arg11b = -(temp3*tempb/arg11) - tempb0 = -(temp4*tempb/temp2) - amixb = tempb0 - tempb1 = -(r*(del1-del2)*temp3*tempb0) - tb = (-(temp1*temp0)-temp4*temp3)*model%r*arvalb + bmix*tempb1 - tempb = arg11b/(del2*b_v+1.0) - b_vb = b_vb + (del1-del2*(del1*b_v+1.0)/(del2*b_v+1.0))*tempb - bmixb = t*tempb1 + b_vb/v - vb = -(bmix*b_vb/v**2) - CALL POPREAL8(bmix) - temp = SUM(n) - nb = nb - bmix*bmixb/temp**2 - bmixb = bmixb/temp - arg10b = 0.0_pr - arg10b = bmixb - arg1b = 0.0_pr - arg1b = amixb - nb = nb + 2*n*model%b*arg10b + 2*n*a**2*arg1b - ab = 0.0_pr - ab = 2*a*n**2*arg1b - CALL POPINTEGER4(ad_to0) - DO i=ad_to0,1,-1 - CALL POPINTEGER4(ad_from) - CALL POPINTEGER4(ad_to) - DO j=ad_to,ad_from,-1 - tempb = (1-kij(i, j))*2*amixb - nij = n(i)*n(j) - nijb = (1-lij(i, j))*(b(i)+b(j))*bmixb + a(i)*a(j)*tempb - ab(i) = ab(i) + nij*a(j)*tempb - ab(j) = ab(j) + nij*a(i)*tempb - nb(i) = nb(i) + n(j)*nijb - nb(j) = nb(j) + n(i)*nijb - END DO - END DO - arg2b = 0.0_pr - WHERE (arg2 .EQ. 0.0) - arg2b = 0.0_pr - ELSEWHERE - arg2b = ab/(2.0*SQRT(arg2)) - END WHERE - result1b = 0.0_pr - result1b = -(k*2*(k*(1.0-result1)+1.0)*ac*arg2b) - arg1(:) = t/tc - arg1b = 0.0_pr - WHERE (arg1 .EQ. 0.0) - arg1b = 0.0_pr - ELSEWHERE - arg1b = result1b/(2.0*SQRT(arg1)) - END WHERE - tb = tb + SUM(arg1b/tc) - arvalb = 0.0_pr - END SUBROUTINE AR_B - - SUBROUTINE AR(model, n, v, t, arval) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:), v, t - REAL(pr), INTENT(OUT) :: arval - REAL(pr) :: amix, a(SIZE(n)), ai(SIZE(n)), z2(SIZE(n)), nij - REAL(pr) :: bmix - REAL(pr) :: b_v - REAL(pr) :: aij(SIZE(n), SIZE(n)), bij(SIZE(n), SIZE(n)) - INTEGER :: i, j - REAL(pr) :: tc(SIZE(n)), pc(SIZE(n)), w(SIZE(n)) - REAL(pr) :: ac(SIZE(n)), b(SIZE(n)), k(SIZE(n)), r - REAL(pr) :: del1, del2 - REAL(pr) :: kij(SIZE(n), SIZE(n)), lij(SIZE(n), SIZE(n)) - INTRINSIC SQRT - INTRINSIC SUM - INTRINSIC LOG - INTRINSIC SIZE - REAL(pr), DIMENSION(SIZE(n)) :: arg1 - REAL(pr), DIMENSION(SIZE(n)) :: result1 - REAL(pr), DIMENSION(SIZE(n)) :: arg2 - REAL(pr), DIMENSION(SIZE(n, 1)) :: arg10 - REAL(pr) :: arg11 - tc = model%tc - pc = model%pc - w = model%w - ac = model%ac - b = model%b - k = model%k - r = model%r - kij = model%kij - lij = model%lij - arg1(:) = t/tc - result1 = SQRT(arg1(:)) - arg2(:) = ac*(1.0+k*(1.0-result1))**2 - a = SQRT(arg2(:)) - amix = 0.0 - bmix = 0.0 - DO i=1,SIZE(n)-1 - DO j=i+1,SIZE(n) - nij = n(i)*n(j) - amix = amix + 2*nij*(a(i)*a(j))*(1-kij(i, j)) - bmix = bmix + nij*(b(i)+b(j))*(1-lij(i, j)) - END DO - END DO - arg1(:) = n**2*a**2 - amix = amix + SUM(arg1(:)) - arg10(:) = n**2*model%b - bmix = bmix + SUM(arg10(:)) - bmix = bmix/SUM(n) - b_v = bmix/v - arg11 = (1.0+del1*b_v)/(1.0+del2*b_v) - arval = (-(SUM(n)*LOG(1.0-b_v))-amix/(r*t*bmix)*1.0/(del1-del2)*LOG(& -& arg11))*(model%r*t) - END SUBROUTINE AR - - PURE FUNCTION VOLUME_INITALIZER(model, n, p, t) RESULT (v0) - IMPLICIT NONE - class(TPR76), INTENT(IN) :: model - REAL(pr), INTENT(IN) :: n(:) - REAL(pr), INTENT(IN) :: p - REAL(pr), INTENT(IN) :: t - REAL(pr) :: v0 - INTRINSIC SUM - v0 = SUM(n*model%b)/SUM(model%b) - END FUNCTION VOLUME_INITALIZER - -END MODULE TAPENADE_PR - diff --git a/fpm.toml b/fpm.toml index 63bad66d4..cd030b676 100644 --- a/fpm.toml +++ b/fpm.toml @@ -38,41 +38,41 @@ main = "demo.f90" [[example]] name = "basics" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "1_basics.f90" [[example]] name = "saturation_points" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "2_saturation_points.f90" [[example]] name = "phase_split" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "3_phase_split.f90" [[example]] name = "phase_envelope" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "4_phase_envelope.f90" [[example]] name = "new_alpha_function" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "new_alpha_function.f90" [[example]] name = "huron_vidal" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "huron_vidal.f90" [[example]] name = "pure_psat" -source-dir = "example/extra" -main = "pure_psat.f90" +source-dir = "example/basics" +main = "5_pure_psat.f90" [[example]] name = "cubic_eos" -source-dir = "example/tutorials" +source-dir = "example/basics" main = "cubic_eos.f90" From d9703998d02b80518d0acbb3e0d7717db05b2906 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 14:08:22 -0300 Subject: [PATCH 10/15] readme to examples --- example/README.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 example/README.md diff --git a/example/README.md b/example/README.md new file mode 100644 index 000000000..f4c3362d8 --- /dev/null +++ b/example/README.md @@ -0,0 +1,14 @@ +# yaeos examples +Here are located a set of examples of how the library could be used. + +We show them in three different categories: + +- basics: The basics of what can be done with the library + - Thermodynamic properties + - Saturation points + - Phase split (FlashPT and FlashVT) + - Phase envelopes tracing + - Implementation of a simple algorithm to calculate pure component + saturation pressure. +- Advanced: How a more advanced user can take benefit of the library. +- extra: Other set of examples with no particular structure. \ No newline at end of file From af4d1c80b9a078056990faec6ea489928ff220a1 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 14:08:30 -0300 Subject: [PATCH 11/15] readme to examples --- example/README.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/example/README.md b/example/README.md index f4c3362d8..9dbf3aeb5 100644 --- a/example/README.md +++ b/example/README.md @@ -4,11 +4,11 @@ Here are located a set of examples of how the library could be used. We show them in three different categories: - basics: The basics of what can be done with the library - - Thermodynamic properties - - Saturation points - - Phase split (FlashPT and FlashVT) - - Phase envelopes tracing - - Implementation of a simple algorithm to calculate pure component + 1. Thermodynamic properties + 2. Saturation points + 3. Phase split (FlashPT and FlashVT) + 4. Phase envelopes tracing + 5. Implementation of a simple algorithm to calculate pure component saturation pressure. - Advanced: How a more advanced user can take benefit of the library. - extra: Other set of examples with no particular structure. \ No newline at end of file From 555ed3101bca0ebf5510776836aec4c9151c9249 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 14:16:41 -0300 Subject: [PATCH 12/15] move some examples --- example/{basics => advanced}/cubic_eos.f90 | 0 example/{basics => advanced}/new_alpha_function.f90 | 0 example/advanced/new_model_adiff.f90 | 0 .../basics/{huron_vidal.f90 => 7_huron_vidal.f90} | 0 example/basics/example_tools.f90 | 5 +---- fpm.toml | 12 ++++++------ 6 files changed, 7 insertions(+), 10 deletions(-) rename example/{basics => advanced}/cubic_eos.f90 (100%) rename example/{basics => advanced}/new_alpha_function.f90 (100%) create mode 100644 example/advanced/new_model_adiff.f90 rename example/basics/{huron_vidal.f90 => 7_huron_vidal.f90} (100%) diff --git a/example/basics/cubic_eos.f90 b/example/advanced/cubic_eos.f90 similarity index 100% rename from example/basics/cubic_eos.f90 rename to example/advanced/cubic_eos.f90 diff --git a/example/basics/new_alpha_function.f90 b/example/advanced/new_alpha_function.f90 similarity index 100% rename from example/basics/new_alpha_function.f90 rename to example/advanced/new_alpha_function.f90 diff --git a/example/advanced/new_model_adiff.f90 b/example/advanced/new_model_adiff.f90 new file mode 100644 index 000000000..e69de29bb diff --git a/example/basics/huron_vidal.f90 b/example/basics/7_huron_vidal.f90 similarity index 100% rename from example/basics/huron_vidal.f90 rename to example/basics/7_huron_vidal.f90 diff --git a/example/basics/example_tools.f90 b/example/basics/example_tools.f90 index 88338a59c..0c9520db5 100644 --- a/example/basics/example_tools.f90 +++ b/example/basics/example_tools.f90 @@ -1,4 +1,4 @@ -! Here there are some functions defined to simplify some examples +!! Here there are some functions defined to simplify some examples module yaeos__example_tools use yaeos @@ -9,13 +9,10 @@ type(CubicEoS) function methane_butane_pr76() result(model) integer, parameter :: nc = 2 real(pr) :: tc(nc), pc(nc), w(nc), kij(nc,nc) - integer :: i - ! Methane/ Butane mixture tc = [190.564, 425.12] ! Critical temperatures pc = [45.99, 37.96] ! Critical pressures w = [0.0115478, 0.200164] ! Acentric factors - model = PengRobinson76(tc, pc, w) end function methane_butane_pr76 end module yaeos__example_tools diff --git a/fpm.toml b/fpm.toml index cd030b676..df2bb31b9 100644 --- a/fpm.toml +++ b/fpm.toml @@ -58,21 +58,21 @@ main = "4_phase_envelope.f90" [[example]] name = "new_alpha_function" -source-dir = "example/basics" +source-dir = "example/advanced" main = "new_alpha_function.f90" [[example]] -name = "huron_vidal" +name = "pure_psat" source-dir = "example/basics" -main = "huron_vidal.f90" +main = "5_pure_psat.f90" [[example]] -name = "pure_psat" +name = "huron_vidal" source-dir = "example/basics" -main = "5_pure_psat.f90" +main = "7_huron_vidal.f90" [[example]] name = "cubic_eos" -source-dir = "example/basics" +source-dir = "example/advanced" main = "cubic_eos.f90" From cff7d2144bd1a9c3dbe01adfba2ba5a6142fbd39 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 15:52:19 -0300 Subject: [PATCH 13/15] disable fixture --- example/advanced/new_alpha_function.f90 | 10 +++++++--- example/basics/2_saturation_points.f90 | 12 ++++++++---- example/basics/4_phase_envelope.f90 | 10 +++++++--- example/basics/example_tools.f90 | 18 ------------------ 4 files changed, 22 insertions(+), 28 deletions(-) delete mode 100644 example/basics/example_tools.f90 diff --git a/example/advanced/new_alpha_function.f90 b/example/advanced/new_alpha_function.f90 index ec261f007..5dd62419e 100644 --- a/example/advanced/new_alpha_function.f90 +++ b/example/advanced/new_alpha_function.f90 @@ -51,11 +51,15 @@ program new_alpha_example type(QMR) :: mr type(MathiasCopeman) :: alpha - real(pr) :: n(2), v, t, P - integer :: i + real(pr) :: n(2), v, t, P, tc(2), pc(2), w(2) + + ! Methane/ Butane mixture + tc = [190.564, 425.12] ! Critical temperatures + pc = [45.99, 37.96] ! Critical pressures + w = [0.0115478, 0.200164] ! Acentric factors ! Get the example PR76 binary model - eos = methane_butane_pr76() + eos = PengRobinson76(tc, pc, w) ! Define the new alpha function parameters alpha%c1 = [0.49258, 0.84209] diff --git a/example/basics/2_saturation_points.f90 b/example/basics/2_saturation_points.f90 index 391c4fa03..92b2e63b1 100644 --- a/example/basics/2_saturation_points.f90 +++ b/example/basics/2_saturation_points.f90 @@ -7,16 +7,20 @@ ! All the outputs of this functions use the `EquilibriumState` type program saturation use yaeos - use yaeos__example_tools, only: methane_butane_pr76 class(ArModel), allocatable :: model type(EquilibriumState) :: sat_point real(pr) :: n(2), T + real(pr) :: Tc(2), Pc(2), w(2) - ! This function returns a predefined model to simplify the step. How it - ! works can be seen in the `example_tools.f90` file. - model = methane_butane_pr76() + ! Methane/ Butane mixture + tc = [190.564, 425.12] ! Critical temperatures + pc = [45.99, 37.96] ! Critical pressures + w = [0.0115478, 0.200164] ! Acentric factors + + ! Get the example PR76 binary model + model = PengRobinson76(tc, pc, w) n = [2.5, 6.7] diff --git a/example/basics/4_phase_envelope.f90 b/example/basics/4_phase_envelope.f90 index 381eb2aaa..725b9e79a 100644 --- a/example/basics/4_phase_envelope.f90 +++ b/example/basics/4_phase_envelope.f90 @@ -8,10 +8,14 @@ program phase_envelope type(PTEnvel2) :: envelope integer, parameter :: nc=2 - real(pr) :: z(nc) + real(pr) :: z(nc), tc(nc), pc(nc), w(nc) - ! Methane/ Butane mixture - model = methane_butane_pr76() + tc = [190.564, 425.12] ! Critical temperatures + pc = [45.99, 37.96] ! Critical pressures + w = [0.0115478, 0.200164] ! Acentric factors + + ! Get the example PR76 binary model + model = PengRobinson76(tc, pc, w) ! Composition z = [0.1_pr, 0.9_pr] diff --git a/example/basics/example_tools.f90 b/example/basics/example_tools.f90 deleted file mode 100644 index 0c9520db5..000000000 --- a/example/basics/example_tools.f90 +++ /dev/null @@ -1,18 +0,0 @@ -!! Here there are some functions defined to simplify some examples -module yaeos__example_tools - use yaeos - -contains - - type(CubicEoS) function methane_butane_pr76() result(model) - !! Binary mixture of methane/butane with PR76 - integer, parameter :: nc = 2 - real(pr) :: tc(nc), pc(nc), w(nc), kij(nc,nc) - - ! Methane/ Butane mixture - tc = [190.564, 425.12] ! Critical temperatures - pc = [45.99, 37.96] ! Critical pressures - w = [0.0115478, 0.200164] ! Acentric factors - model = PengRobinson76(tc, pc, w) - end function methane_butane_pr76 -end module yaeos__example_tools From 7a9ec6fcde2c34390081b3b152c1b56c1811bf46 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 15:58:21 -0300 Subject: [PATCH 14/15] fix --- example/basics/4_phase_envelope.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/example/basics/4_phase_envelope.f90 b/example/basics/4_phase_envelope.f90 index 725b9e79a..2e1be8d19 100644 --- a/example/basics/4_phase_envelope.f90 +++ b/example/basics/4_phase_envelope.f90 @@ -1,6 +1,5 @@ program phase_envelope use yaeos - use yaeos__example_tools, only: methane_butane_pr76 implicit none class(ArModel), allocatable :: model ! Model From 96a3e9a97ae2863b7806aaab03a679f6d550f780 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Fri, 19 Jul 2024 16:00:58 -0300 Subject: [PATCH 15/15] fix --- example/advanced/new_alpha_function.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/example/advanced/new_alpha_function.f90 b/example/advanced/new_alpha_function.f90 index 5dd62419e..529613cf7 100644 --- a/example/advanced/new_alpha_function.f90 +++ b/example/advanced/new_alpha_function.f90 @@ -44,7 +44,6 @@ end subroutine alpha end module alpha_mathias_copeman program new_alpha_example - use yaeos__example_tools, only: methane_butane_pr76 use yaeos, only: pr, PengRobinson76, CubicEoS, QMR use alpha_mathias_copeman, only: MathiasCopeman type(CubicEoS) :: eos