possibility to only calculate fugacity coef

ipqa-research · Feb 5, 2024 · 6ef685f · 6ef685f
1 parent c69cfac
commit 6ef685f
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 6 deletions.
diff --git a/example/benchmark.f90 b/example/benchmark.f90
@@ -46,7 +46,7 @@ subroutine yaeos_run(n, dn, model_name)
         if (dn) then
             call fugacity_vt(model, z, V, T, P, lnfug, dlnPhidP, dlnphidT, dlnphidn)
         else
-            call fugacity_vt(model, z, V, T, P, lnfug, dlnPhidP, dlnphidT)
+            call fugacity_vt(model, z, V, T, lnfug=lnfug)
         end if
     end subroutine
 
@@ -57,13 +57,13 @@ subroutine benchmarks
         real(8) :: et, st
 
 
-        do n=1,5
+        do n=1,30
             time = 0
             std = 0
             mean = 0
             do i=1,nevals
                 call cpu_time(st)
-                    call yaeos_run(n, .true., "Adiff PR76")
+                    call yaeos_run(n, .true., "Analytic PR76")
                 call cpu_time(et)
 
                 time = (et-st)*1e6

diff --git a/src/thermoprops.f90 b/src/thermoprops.f90
@@ -62,7 +62,16 @@ subroutine fugacity_vt(self, &
       RT = R*T
       Z = V/(TOTN*RT) ! this is Z/P
 
-      if (present(dlnphidn)) then
+      if (present(lnfug) .and. .not. (&
+         present(dlnphidn) &
+         .and. present(dlnphidp) &
+         .and. present(dlnphidt) &
+         .and. present(p) &
+      )) then
+         call self%residual_helmholtz(n, v, t, Arn=Arn)
+         lnfug(:) = Arn(:)/RT - log(Z)
+         return
+      else if (present(dlnphidn)) then
          call self%residual_helmholtz(&
             n, V, T, Ar=Ar, ArV=ArV, ArV2=ArV2, ArTV=ArTV, &
             Arn=Arn, ArVn=ArVn, ArTn=ArTn, Arn2=Arn2 &
@@ -73,14 +82,14 @@ subroutine fugacity_vt(self, &
             Arn=Arn, ArVn=ArVn, ArTn=ArTn &
          )
       end if
+
+      lnfug(:) = Arn(:)/RT - log(Z)
 
       P = TOTN*RT/V - ArV
       dPdV = -ArV2 - RT*TOTN/V**2
       dPdT = -ArTV + TOTN*R/V
       dPdN(:) = RT/V - ArVn(:)
 
-      lnfug(:) = Arn(:)/RT - log(Z)
-
       if (present(dlnphidp)) dlnphidp(:) = -dPdN(:)/dPdV/RT - 1._pr/P
       if (present(dlnphidt)) then
          dlnphidt(:) = (ArTn(:) - Arn(:)/T)/RT + dPdN(:)*dPdT/dPdV/RT + 1._pr/T