diff --git a/.gitignore b/.gitignore index 32b8afb..a01e08b 100644 --- a/.gitignore +++ b/.gitignore @@ -86,4 +86,14 @@ infiles/YPF10 YPF10 fig77 fpm.toml - +src/tools/ftools__io.mod +YPF1.nml +fpm.toml +src/molecular_weigth/molecular_weigth.mod +.gitignore +src/bfr_routines/bfr_routines.mod +src/characterization/characterization.mod +src/critical_parameters/critical_parameters.mod +.gitignore +src/density/density.mod +src/lumping/lumping.mod diff --git a/PVT1.nml b/PVT1.nml new file mode 100644 index 0000000..5340f70 --- /dev/null +++ b/PVT1.nml @@ -0,0 +1,72 @@ +! Namelist based input file +! ============================================================================== +! - PVT1 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0088 0.0052 0.6037 0.1259 0.0728 0.0146 0.0308 0.0117 0.0143 + ! compositions of defined components + scn_z = 0.0189 0.0154 0.0192 0.0121 0.0088 0.0067 0.0054 0.0048 0.0037 0.0029 0.0022 0.0018 0.0016 0.0013 0.0010 0.0008 0.0007 0.0006 0.0005 0.0004 0.0004 0.0003 0.0003 0.0002 + ! compositions of single carbon numbers + plus_z = 0.0022 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1122 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.0933 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0311 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0022 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 86.2 96.4 108.7 123.9 138.3 146.9 160.3 174.8 190.7 207.8 223.7 241.6 247.9 261.1 267 290.8 299.7 298.8 320.4 337.4 350.7 364 377.3 390.6 + ! molecular weights of single carbon numbers + plus_mw = 561.11 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 146.72 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 157.58 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 243.55 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 561.11 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9816 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.7969 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8113 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8596 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9816 + ! density of residual fraction CN 30 plus +/ diff --git a/PVT2.nml b/PVT2.nml new file mode 100644 index 0000000..42be41b --- /dev/null +++ b/PVT2.nml @@ -0,0 +1,73 @@ +! Namelist based input file +! ============================================================================== +! - PVT2 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0021 0.0035 0.6048 0.1317 0.0604 0.0125 0.0279 0.0120 0.0148 + ! compositions of defined components + scn_z = 0.0188 0.0186 0.0233 0.0149 0.0107 0.0076 0.0057 0.0054 0.0041 0.0030 0.0024 0.0020 0.0018 0.0016 0.0012 0.0011 0.0009 0.0008 0.0007 0.0006 0.0006 0.0005 0.0004 0.0004 + ! compositions of single carbon numbers + plus_z = 0.0032 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1303 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.1114 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0364 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0032 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 85.92 95.74 107.97 123.86 138.88 146.65 162.24 174.24 189.11 207.29 220.29 238.25 255.75 262.48 282.67 279 309.59 312.96 334.59 346.06 359.99 373.92 387.86 401.79 + ! molecular weights of single carbon numbers + plus_mw = 553.6 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 150.4 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 159.93 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 253.63 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 553.6 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9807 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.8022 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8146 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8654 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9807 + ! density of residual fraction CN 30 plus +/ + diff --git a/PVT4.nml b/PVT4.nml new file mode 100644 index 0000000..d3a0602 --- /dev/null +++ b/PVT4.nml @@ -0,0 +1,72 @@ +! Namelist based input file +! ============================================================================== +! - PVT4 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0054 0.0023 0.7261 0.1179 0.0532 0.0112 0.0176 0.011 0.0074 + ! compositions of defined components + scn_z = 0.00978 0.00778 0.01038 0.00559 0.00369 0.00249 0.0017 0.0016 0.0011 0.0009 0.0006 0.0005 0.0004 0.0003 0.00021 0.00016 0.00013 0.0001 0.00008 0.00007 0.00005 0.00004 0.00003 0.00003 + ! compositions of single carbon numbers + plus_z = 0.0002 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.0478 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.0381 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0082 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0002 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 85.71 96.64 108.98 124.03 137.52 145.38 161.53 169.94 186.01 197.43 233.32 242.3 255.76 269.22 279.99 294.34 308.69 323.04 337.39 351.74 366.08 380.43 394.78 409.13 + ! molecular weights of single carbon numbers + plus_mw = 551.2 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 126.64 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 134.38 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 217.35 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 551.2 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9804 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.7728 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.7889 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8414 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9804 + ! density of residual fraction CN 30 plus +/ diff --git a/PVT5.nml b/PVT5.nml new file mode 100644 index 0000000..dbcaad7 --- /dev/null +++ b/PVT5.nml @@ -0,0 +1,74 @@ +! Namelist based input file +! ============================================================================== +! - PVT5 from whitson dataset +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 24 + ! number of single carbon number components + scn_nc_ps = 30 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 + ! names of scn fractions + scn_plus= "C30+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.0026 0.0053 0.5403 0.1260 0.0764 0.0139 0.0341 0.0136 0.0168 + ! compositions of defined components + scn_z = 0.0213 0.0183 0.0205 0.0152 0.0112 0.0098 0.0079 0.0072 0.0055 0.0040 0.0031 0.0027 0.0025 0.0022 0.0017 0.0015 0.0013 0.0012 0.0010 0.0010 0.0008 0.0007 0.0006 0.0006 + ! compositions of single carbon numbers + plus_z = 0.0057 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.1711 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.1492 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.0512 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.0057 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.2 72.2 + ! molecular weights of defined components + scn_mw = 86 95.9 108.2 124 138.3 147.5 161.1 174.1 189.5 207.7 223.1 236.6 247.8 262 276.8 297.8 310.3 312.2 325.3 353.6 360.2 377.3 384.0 401.5 + ! molecular weights of single carbon numbers + plus_mw = 579.2 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 157.21 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 167.73 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 265.89 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 579.2 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + ! density of single carbon numbers + plus_density = 0.9837 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.8102 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.8214 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.8732 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.9837 + ! density of residual fraction CN 30 plus +/ + + diff --git a/YPF2.nml b/YPF2.nml new file mode 100644 index 0000000..d8ffef1 --- /dev/null +++ b/YPF2.nml @@ -0,0 +1,78 @@ +! Namelist based input file +! ============================================================================== +! - YPF2 +! - Composition: molar +! - Molecular Weights: g/mol +! - Density: g/cm3 +! +! +&nml_setup + def_comp_nc = 9 + ! number of defined components before scn + scn_nc = 14 + ! number of single carbon number components + scn_nc_ps = 20 + ! CN from which all SCN fractions will be lumped into the specified number of pseudos + numbers_ps = 5 + ! number of pseudos in which the scn fractions grouped + density_setup = .true. + ! Logical variable to choose the density calculation method + number_plus_density = 6 + ! CN from which we considered ....... +/ +&nml_components + def_components = "N2" "CO2" "C1" "C2" "C3" "iC4" "nC4" "iC5" "nC5" + ! names of defined components before scn + scn = 6 7 8 9 10 11 12 13 14 15 16 17 18 19 + ! names of scn fractions + scn_plus= "C20+" + ! name of plus fraction +/ +&nml_composition + def_comp_z = 0.00005 0.01963 0.06271 0.00707 0.01354 0.0045 0.01607 0.01069 0.01562 + ! compositions of defined components + scn_z = 0.03374 0.06051 0.08007 0.06365 0.05495 0.03993 0.03139 0.03392 0.03133 0.03183 0.03042 0.02984 0.03127 0.02964 + ! compositions of single carbon numbers + plus_z = 0.26762 + ! composition of residual fraction CN 20 plus + plus6_z_exp = 0.26762 + ! composition of residual fraction CN 6 plus + plus7_z_exp = 0.26762 + ! composition of residual fraction CN 7 plus + plus12_z_exp = 0.26762 + ! composition of residual fraction CN 12 plus + plus30_z_exp = 0.26762 + ! composition of residual fraction CN 30 plus +/ +&nml_molecular_weight + def_comp_mw = 28.01 44.01 16 30.1 44.1 58.1 58.1 72.1 72.1 + ! molecular weights of defined components + scn_mw = 84 96 107 121 134 147 161 175 190 206 222 237 251 263 + ! molecular weights of single carbon numbers + plus_mw = 543.35 + ! molecular weight of residual fraction CN 20 plus + plus6_mw_exp = 465 + ! molecular weight of residual fraction CN 6 plus + plus7_mw_exp = 500 + ! molecular weight of residual fraction CN 7 plus + plus12_mw_exp = 600 + ! molecular weight of residual fraction CN 12 plus + plus30_mw_exp = 700 + ! molecular weight of residual fraction CN 30 plus +/ +&nml_density + scn_density = 0.685 0.749 0.768 0.793 0.808 0.815 0.836 0.85 0.861 0.873 0.882 0.873 0.875 0.885 + ! density of single carbon numbers + plus_density = 0.936 + ! density of residual fraction CN 20 plus + plus6_density_exp = 0.77 + ! density of residual fraction CN 6 plus + plus7_density_exp = 0.77 + ! density of residual fraction CN 7 plus + plus12_density_exp = 0.77 + ! density of residual fraction CN 12 plus + plus30_density_exp = 0.77 + ! density of residual fraction CN 30 plus +/ + + diff --git a/app/main.f90 b/app/main.f90 index 170cbc0..bc7f772 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -1,21 +1,23 @@ program main - use data_from_input, only: data_from_file, FluidData - use routines + !use data_from_input, only: data_from_file, FluidData + use characterization implicit none - integer :: k, i + type(FluidDataOut) :: prueba - type(FluidData) :: fluid + integer :: i + + prueba = characterize(file='PVT5.nml', mw_source="experimental", method = "plus_mw", pho_method = 3 , fix_C=.true., eos='PR') + write(*, *) prueba + print*, prueba%C, prueba%a, prueba%b, prueba%plus_mw, prueba%plus_z + print*, prueba%c_max + !print*, prueba%a , prueba%b , prueba%n_init, prueba%input_data%number_plus_density + !do i = 1, size(prueba%mol_fraction) + ! print*, prueba%mol_fraction(i) + !end do - prueba = characterize(file='oil1.nml', mw_source="calculated", method = "plus_mw", fix_C=.true.) - print*, prueba%n_init - print*, prueba%a , prueba%b - print*, prueba%C, prueba%plus_mw - print*, prueba%a_d, prueba%b_d - print*, "-------------------------------------------------------------------" - print*, prueba%nc_plus end program main diff --git a/app/optimize.f90 b/app/optimize.f90 new file mode 100644 index 0000000..ceeaacc --- /dev/null +++ b/app/optimize.f90 @@ -0,0 +1,65 @@ +module my_objective + use ForTimize, only: pr + + type :: ExpData + real(pr), allocatable :: x(:) + real(pr), allocatable :: y(:) + end type + +contains + + subroutine foo(x, F, dF, data) + real(pr), intent(in) :: x(:) + !! Vector of parameters to optimize + real(pr), intent(out) :: F + !! Value of the objetive function after optimization. + real(pr), optional, intent(out) :: dF(:) !! + !! Optional gradient. + class(*), optional, intent(in out) :: data + !! Special data that the function could use. + + real(pr) :: a, b, c + + a = x(1) + b = x(2) + c = x(3) + + F = 0 + + select type(data) + type is (ExpData) + F = sum( (data%y - (quad(data%x, a, b, c)))**2 ) + print *, F, data%x + end select + end subroutine + + elemental function quad(x, a, b, c) + real(pr), intent(in) :: x, a, b, c + real(pr) :: quad + + quad = a*x**2 + b*x + c + end function +end module + + +program oscar + use ForTimize, only: pr, minimize + use my_objective, only: foo, ExpData + + real(pr) :: x(3), F + type(ExpData) :: exp + + exp%x = [2, 5, 7] + exp%y = [4, 22, 50] + + ! Initial guess + x = [1, 4, 3] + + ! Minimize uses the Nelder-Mead algorithm as a default + call minimize(foo, x, F, data=exp) + ! call foo(x, F, data=exp) + + ! Print results + print *, x + print *, F +end program oscar diff --git a/app/pedersen.f90 b/app/pedersen.f90 new file mode 100644 index 0000000..3975ebb --- /dev/null +++ b/app/pedersen.f90 @@ -0,0 +1,76 @@ +module my_objective1 + use ForTimize, only: pr + + + type :: ExpData + real(pr), allocatable :: x(:) + real(pr), allocatable :: zplus(:) + real(pr), allocatable :: mplus(:) + end type ExpData + +contains + + subroutine foo(x, F, dF, data) + real(pr), intent(in) :: x(:) + !! Vector of parameters to optimize + real(pr), intent(out) :: F + !! Value of the objetive function after optimization. + real(pr), optional, intent(out) :: dF(:) !! + !! Optional gradient. + class(*), optional, intent(in out) :: data + !! Special data that the function could use. + + real(pr) :: a, b + integer :: i + integer, parameter :: n = 181 + real(pr) :: cn(n) + real(pr) :: z_i(n) + real(pr) :: m_i(n) + real(pr) :: zplus_cal + real(pr) :: mplus_cal + + + do i = 1, n + cn(i) = 20.0d0 + i - 1 + end do + + a = x(1) + b = x(2) + + z_i = exp(a*cn+b) + m_i = 14*cn-4 + zplus_cal = sum(z_i) + mplus_cal = sum(z_i*m_i)/zplus_cal + + F = 0 + + select type(data) + type is (ExpData) + F = sum((data%zplus - (zplus_cal))**2 + (data%mplus - (mplus_cal))**2) + end select + end subroutine foo + +end module my_objective1 + + +program oscar + use ForTimize, only: pr, minimize + use my_objective1, only: foo, ExpData + + real(pr) :: x(2), F + type(ExpData) :: exp + + exp%zplus = [0.1107] + exp%mplus = [419.3] + + ! Initial guess + x = [ -9.3214408744470914E-002, -2.7557414112453325 ] + + ! Minimize uses the Nelder-Mead algorithm as a default + call minimize(foo, x, F, data=exp) + + ! Print results + print *, x + print *, F + !print*, cn +end program oscar diff --git a/fpm.toml b/fpm.toml index 77a72ab..d4d70ef 100644 --- a/fpm.toml +++ b/fpm.toml @@ -14,4 +14,6 @@ library = false [fortran] implicit-typing = false implicit-external = false -source-form = "free" \ No newline at end of file +source-form = "free" +[dependencies] +ForTimize = {git='https://github.com/fedebenelli/ForTimize'} \ No newline at end of file diff --git a/oil1.nml b/oil1.nml index 081200a..61143a2 100644 --- a/oil1.nml +++ b/oil1.nml @@ -49,3 +49,11 @@ ! density of residual fraction / + +&nml_density + scn_density = 0.685 0.749 0.768 0.793 0.808 0.815 0.836 0.85 0.861 0.873 0.882 0.873 0.875 0.885 + ! density of single carbon numbers + plus_density = 0.936 + ! density of residual fraction +/ + diff --git a/src/bfr_routines/bfr_routines.f90 b/src/bfr_routines/bfr_routines.f90 new file mode 100644 index 0000000..5010f9d --- /dev/null +++ b/src/bfr_routines/bfr_routines.f90 @@ -0,0 +1,226 @@ +module bfr_routines + use constants + +contains + + subroutine Linear_Regression(x, y, a, b, r2) + !! This subroutine computes the regression line for a data set of x, y variables. + implicit none + + integer :: i + integer :: n !! total number of data set. + real(pr), intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable + real(pr), intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable + real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables + real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line + real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line + real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient + + n = size(x) + !Calculation of a y b --> y=a*x+b + t1=0.;t2=0.;t3=0.;t4=0.; + + do i=1,n + t1=t1+x(i)*y(i) + t2=t2+x(i) + t3=t3+y(i) + t4=t4+x(i)**2 + end do + + a=(n*t1-t2*t3)/(n*t4-t2**2) + b=(t3-a*t2)/n + !coefficient calculation of correlations r2 + aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; + + do i=1, n + aux1= aux1 + x(i)*y(i) + aux2= aux2 + x(i) + aux3= aux3 + y(i) + aux4= aux4 + x(i)**2 + aux5= aux5 + y(i)**2 + end do + + r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) + + end subroutine Linear_Regression + + subroutine Best_Linear_Regression(scn_nc, scn, scn_z, plus_z, a, b, r2, & + n_init, c_max_blr) + !! This subroutine calculates the best regression line for an fluid. + implicit none + + integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid + integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid + integer, intent(out) :: c_max_blr !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) + real(pr), intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts + real(pr), intent(in) :: plus_z !! composition of residual fraction from input file + real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. + real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. + real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. + integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression + integer :: i, j, k , k_old, n_best, x_aux ! internal variables + real(pr), dimension(scn_nc) :: x_blr, y_blr ! x_blr: carbon number vector, y_blr: logarithm of mole fraction vector; for each step + real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux ! auxiliary variables + + k=5 + r2=0.0001_pr + r2_old=0.00001_pr + r2_best=0.0001_pr + + do while (r2.gt.r2_old.or.r2_old.lt.0.9) + k_old=k + r2_old=r2 + a_old=a + b_old=b + + if (r2.gt.r2_best) then + r2_best=r2 + a_best=a + b_best=b + n_best=scn(scn_nc-k+2) + end if + + if (k.gt.scn_nc) then + if (r2.gt.r2_best)then + n_init=scn(scn_nc-k+2) + go to 22 + else + r2=r2_best + a=a_best + b=b_best + n_init=n_best + go to 22 + end if + end if + + j=1 + x_blr=0. + y_blr=0. + + do i=scn_nc-k+1, scn_nc + x_blr(j)=scn(i) + y_blr(j)=scn_z(i) + j=j+1 + end do + + call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) + k=k+1 + end do + + r2=r2_old + a=a_old + b=b_old + n_init=scn(scn_nc-k_old+2) + + 22 continue + + z_sum = 0d0 + x_aux=scn(scn_nc) + + do while (z_sum.lt.plus_z.and.x_aux<300) + x_aux=x_aux+1 + z_aux= exp(a*x_aux+b) + z_sum=z_sum+z_aux + end do + + c_max_blr=x_aux + + !print*, a, b, c_max_blr, r2, n_init + + end subroutine Best_Linear_Regression + + subroutine LimitLine(scn_nc, scn, plus_z, a_blr, b_blr, a_lim, b_lim, & + c_max_lim, half) + !! This subroutine obtains the limit line constants for an fluid. + implicit none + + integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid. + integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid. + integer, intent(out) :: c_max_lim !! maximum carbon number obtained for limit line. + !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr). + real(pr), intent(in) :: a_blr !! input constant from the Best linear regression. + real(pr), intent(in) :: b_blr !! input constant from the Best linear regression. + real(pr), intent(in) :: plus_z !! composition of residual fraction from input file. + real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line. + real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. + real(pr), intent(out) :: half ! variable defined in the numerical method to converge the area under the curve of the function. + real(pr) :: z_lim, cross_cn, z_cross, z_aux + integer :: x_aux ! internal variable + + !A and B limit calculation. + z_lim = 0.0_pr + half = 0.506_pr ! modified by oscar, the variable half was removed from the argument + + do while (z_lim.lt.plus_z) + half = half + 0.002d0 + cross_cn = scn(scn_nc)+half ! typically 19.508 in first try + z_cross = exp((a_blr*cross_cn) + b_blr) + a_lim = -(z_cross)/plus_z + b_lim = log(z_cross)-(a_lim*cross_cn) + x_aux=scn(scn_nc) + z_lim = 0._pr + do while (z_lim.lt.plus_z.and.x_aux<300) + x_aux=x_aux+1 + z_aux= exp(a_lim*x_aux+b_lim) + z_lim=z_lim+z_aux + end do + continue + end do + + c_max_lim = x_aux + + end subroutine LimitLine + + subroutine Line_C60_max (scn_nc, scn, plus_z, a_blr, b_blr, half, a_60, & + b_60, c_max_60) + !! This subroutine obtains the Cmax60 line constants for an fluid. + + implicit none + integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid + integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid + integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. + real(pr), intent(in) :: a_blr !! input constant from the Best linear regression + real(pr), intent(in) :: b_blr !! input constant from the Best linear regression + real(pr), intent(in) :: plus_z !! composition of residual fraction from input file + real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line + real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. + real(pr), intent(in) :: half !variable defined in the numerical method to converge the area under the curve of the function + integer :: x_aux ! auxiliary variable + real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range ! internal variables + real(pr) :: a_old, Zp_60, F, dF_dA ! internal variables + + cross_cn = scn(scn_nc)+half ! typically 19.51 + z_cross= exp((a_blr*cross_cn)+b_blr) + F_tol=1_pr + a_tol=1_pr + a_60=a_blr + var_range = 60.5_pr - cross_cn + + do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) + a_old=a_60 + Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 + F = plus_z - Zp_60 + dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected + a_60=a_old-(F/dF_dA) + a_tol=abs(a_old-a_60) + F_tol=abs(F) + end do + + b_60=log(z_cross)-a_60*cross_cn + z_sum=0.0_pr + x_aux=scn(scn_nc) + + do while (z_sum.lt.plus_z) + x_aux=x_aux+1 + z_aux= exp(a_60*x_aux+b_60) + z_sum=z_sum+z_aux + end do + + c_max_60=x_aux + + !print*, a_60,b_60,c_max_60 + end subroutine Line_C60_max + + +end module bfr_routines + diff --git a/src/characterization/characterization.f90 b/src/characterization/characterization.f90 new file mode 100644 index 0000000..9d42ae0 --- /dev/null +++ b/src/characterization/characterization.f90 @@ -0,0 +1,58 @@ +module characterization + + use constants + use dtypes, only: FluidData, FluidDataOut + use data_from_input, only: data_from_file, FluidData + use molecular_weigth + use density + use lumping + use critical_parameters + + +contains + + type(FluidDataOut) function characterize(file, mw_source, method, pho_method, fix_C, eos) & + result(characterization) + + type(FluidData) :: fluid + character(len=*), intent(in) :: file !! file name + character(len=*), intent(in), optional :: method !! plus_mw or global_mw + character(len=*), intent(in) :: mw_source + logical, intent(in), optional :: fix_C + character(len=*), intent(in) :: eos + integer, intent(in) :: pho_method !! selector method of density funtion. + + + + fluid = data_from_file(file) + characterization%input_data = fluid + allocate(characterization%scn_z(fluid%scn_nc)) + allocate(characterization%log_scn_z(fluid%scn_nc)) + allocate(characterization%scn_mw(fluid%scn_nc)) + allocate(characterization%carbon_number_plus(0)) + allocate(characterization%plus_z_i(0)) + allocate(characterization%product_z_mw_plus_i(0)) + allocate(characterization%scn_zm(fluid%scn_nc)) + allocate(characterization%scn_i(0)) + allocate(characterization%plus6_density(0)) + allocate(characterization%mol_fraction(0)) + allocate(characterization%lumped_z(0)) + allocate(characterization%lumped_mw(0)) + allocate(characterization%lumped_densities(0)) + allocate(characterization%critical_temperature(0)) + allocate(characterization%critical_pressure(0)) + allocate(characterization%acentric_factor(0)) + allocate(characterization%m_funtion(0)) + + + call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) + call density_funtion(fluid=fluid, mw_source=mw_source, pho_method= pho_method, characterization=characterization) + call lump(fluid=fluid, characterization=characterization) + call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) + + + end function characterize + +end module characterization + + diff --git a/src/critical_parameters/critical_parameters.f90 b/src/critical_parameters/critical_parameters.f90 new file mode 100644 index 0000000..39bfcf4 --- /dev/null +++ b/src/critical_parameters/critical_parameters.f90 @@ -0,0 +1,72 @@ +module critical_parameters + + use constants + use dtypes + use defined_critical_parameters + +contains + + subroutine get_critical_constants(fluid, characterization, eos) + !! this subroutine doing: + !! - Calculation of Tc, Pc, omega from Pedersen correlations (adapted from code "C7plusPedersenTcPcOm") + !! - Calculation of EoS parameters (taken from CubicParam) + !use critical_parameters + + implicit none + type(FluidData), intent(in) :: fluid + type(FluidDataOut), intent(inout) :: characterization + character(len=*), intent(in) :: eos + real(pr), allocatable :: modified_c(:) + real(pr), allocatable :: tc(:) + real(pr), allocatable :: pc(:) + real(pr), allocatable :: om(:) + real(pr), allocatable :: m_funtion(:) + integer :: i + + associate(& + rho => characterization%lumped_densities, & + mw => characterization%lumped_mw & + ) + + allocate(modified_c(0)) + allocate(tc(0)) + allocate(pc(0)) + allocate(om(0)) + allocate(m_funtion(0)) + + call get_parameteres_for_critical(eos) + + tc = c1*rho + c2*log(mw) + c3*mw + (c4/mw) + pc = exp(d1_p + d2*(rho**(d5)) + (d3/mw) + (d4/(mw**2))) + + if (eos == "SRK") then + m_funtion = e1 + e2*mw + e3*rho + e4*(mw**(2)) + else + if (eos == "PR" .or. eos =="RKPR") then + m_funtion = e3*rho + 0.2 + 1.8*(1-exp(-mw/220)) + ! Alternative expression without max with MW - May 2019 + endif + endif + + modified_c = c0 - m_funtion + om = (-b + sqrt(b**2 - 4*a*modified_c))/(2*a) + + !do i = 1, size(mw) + ! print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) + !end do + + characterization%critical_temperature = tc + characterization%critical_pressure = pc + characterization%acentric_factor = om + characterization%m_funtion = m_funtion + + end associate + + end subroutine get_critical_constants + +end module critical_parameters + + + + + diff --git a/src/critical_parameters/defined_critical_param.f90 b/src/critical_parameters/defined_critical_param.f90 new file mode 100644 index 0000000..b00e1d1 --- /dev/null +++ b/src/critical_parameters/defined_critical_param.f90 @@ -0,0 +1,82 @@ +module defined_critical_parameters + use constants + + implicit none + + real(pr), allocatable :: tc_def(:), pc_def(:), om_def(:) + real(pr):: c1, c2, c3, c4, d1_p, d2, d3, d4, d5, e1, e2, e3, e4 + real(pr):: a, b, c0, del_1 + +contains + + subroutine get_parameteres_for_critical(eos) + implicit none + character(len=*) :: eos + + tc_def = [126.2, 304.21, 190.564, 305.32, 369.83, 408.14, 425.12, 460.43, 469.7] + pc_def = [34.0, 73.83, 45.99, 48.72, 42.48, 36.48, 37.96, 33.81, 33.7] + om_def = [0.038, 0.224, 0.012, 0.099, 0.152, 0.181, 0.20, 0.228, 0.252] + + select case(eos) + + case("SRK") + c1 = 1.6312d2 + c2 = 8.6052d1 + c3 = 4.3475d-1 + c4 = -1.8774d3 + d1_p = -1.3408d-1 + d2 = 2.5019 + d3 = 2.0846d2 + d4 = -3.9872d3 + d5 = 1.0d0 + e1 = 7.4310d-1 + e2 = 4.8122d-3 + e3 = 9.6707d-3 + e4 = -3.7184d-6 + a = -0.176 + b = 1.574 + c0 = 0.48 + !del_1 = 1.0D0 + + case("PR") ! Coeficientes originales Tabla 5.3 Pedersen + c1 = 7.34043d1 + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1_p = 7.28462d-2 + d2 = 2.18811d0 + d3 = 1.63910d2 + d4 = -4.04323d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + a = -0.26992 + b = 1.54226 + c0 = 0.37464 + !del_1 = 1.0D0 + sqrt(2.0) + + case("RKPR") + c1 = 7.34043d1 ! Coeficientes originales Tabla 5.3 Pedersen + c2 = 9.73562d1 + c3 = 6.18744d-1 + c4 = -2.05932d3 + d1_p = 2.3d0 + d2 = -0.5d0 + d3 = 1.85d2 + d4 = -4.0d3 + d5 = 0.25d0 + e1 = 3.73765d-1 + e2 = 5.49269d-3 + e3 = 1.17934d-2 + e4 = -4.93049d-6 + a = -0.26992 + b = 1.54226 + c0 = 0.37464 + + end select + + end subroutine get_parameteres_for_critical + +end module defined_critical_parameters \ No newline at end of file diff --git a/src/density/density_funtion.f90 b/src/density/density_funtion.f90 new file mode 100644 index 0000000..32eb4b1 --- /dev/null +++ b/src/density/density_funtion.f90 @@ -0,0 +1,160 @@ +module density + use constants + use dtypes, only: FluidData, FluidDataOut + +contains + + subroutine density_funtion(fluid, mw_source, pho_method, characterization) + !! this subroutine ... + implicit none + type(FluidData) :: fluid + type(FluidDataOut) :: characterization + character(len=*), intent(in) :: mw_source + integer, intent(in) :: pho_method + !real(pr), allocatable :: volume_6plus_cal(:) + real(pr) :: volume_exp + real(pr) :: a_d_old, aux + real(pr) :: difference, difference_old + + + !volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & + ! (fluid%product_z_mw_plus/fluid%plus_density) + + + + + + + associate( & + a_d => characterization%a_d , & + b_d => characterization%b_d, & + volume_cal => characterization%volume_cal, & + volume_exp => characterization%volume_exp & + ) + + call density_method (fluid, pho_method,characterization) + + ! Now find constants for the Density function + a_d = -0.50_pr ! initial guess + b_d = 0.685_pr - a_d*exp(-0.60_pr) + call calculate_volume(fluid, mw_source, pho_method, characterization) + difference_old = volume_cal - volume_exp + a_d_old = a_d + a_d = -0.49_pr + b_d = 0.685_pr - a_d*exp(-0.60_pr) + call calculate_volume(fluid, mw_source, pho_method, characterization) + difference = volume_cal - volume_exp + + do while (abs(difference) > 0.001_pr) + aux = a_d + a_d = a_d - difference*(a_d - a_d_old)/(difference - difference_old) + b_d = 0.685_pr - a_d*exp(-0.60_pr) + a_d_old = aux + difference_old = difference + call calculate_volume(fluid, mw_source, pho_method, characterization) + difference = volume_cal - volume_exp + end do + + end associate + end subroutine density_funtion + + subroutine density_method (fluid, pho_method,characterization) + !! this subroutine ... + implicit none + type(FluidData) :: fluid + type(FluidDataOut) :: characterization + integer, intent(in) :: pho_method + real(pr) :: volume_exp + + select case (pho_method) ! cambiar a rho + + case (1) + ! case 1 use experimental scn's densities reported to calculate + ! experimental volume for 6plus. The densities values are calculated + ! since CN plus + volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & + (fluid%product_z_mw_plus/fluid%plus_density) + characterization%volume_exp = volume_exp + + case (2) + !This case use scn densities reported to calculate experimental volume + ! this densities doesn't experimental, they can be reference values. + ! the densities values are calculated since 6 plus + + volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & + (fluid%product_z_mw_plus/fluid%plus_density) + characterization%volume_exp = volume_exp + + case (3) + ! This case use one experimental value (6plus) to calculate the density + ! distribution + + volume_exp = (fluid%plus6_z_exp*fluid%plus6_mw_exp)/(fluid%plus6_density_exp) + characterization%volume_exp = volume_exp + + end select + + + + + + + + + end subroutine density_method + + subroutine calculate_volume(fluid, mw_source, pho_method, characterization) + !! this subroutine ... + implicit none + type(FluidData), intent(in) :: fluid + type(FluidDataOut), intent(inout) :: characterization + character(len=*), intent(in) :: mw_source + integer, intent(in) :: pho_method + real(pr) :: volume_cal + real(pr), allocatable :: plus_density_cal(:), scn_density_cal(:), plus6_density(:) + + !allocate(density_plus_cal(0), density_scn_cal(0), plus6_density(0)) + + + associate(& + a_d => characterization%a_d, & + b_d => characterization%b_d, & + cn_plus => characterization%carbon_number_plus, & + z_m_plus_i => characterization%product_z_mw_plus_i, & + scn_z_i => characterization%scn_z, & + scn_mw_i => characterization%scn_mw, & + scn_zm => characterization%scn_zm & + ) + + scn_density_cal = ((a_d)*(exp(- (real(fluid%scn, pr))/10._pr))) + b_d + plus_density_cal =((a_d)*(exp(- (real(cn_plus, pr) )/10._pr))) + b_d + + + select case (pho_method) + + case(1) + volume_cal = sum((scn_zm)/(fluid%scn_density)) + & + sum(z_m_plus_i/plus_density_cal) + plus6_density = [fluid%scn_density, plus_density_cal] + + case(2) + volume_cal = sum((scn_zm)/(scn_density_cal)) + & + sum(z_m_plus_i/plus_density_cal) + plus6_density = [scn_density_cal, plus_density_cal] + + case(3) + volume_cal = sum((scn_zm)/(scn_density_cal)) + & + sum(z_m_plus_i/plus_density_cal) + plus6_density = [scn_density_cal, plus_density_cal] + + end select + + characterization%volume_cal = volume_cal + characterization%plus6_density = plus6_density + + end associate + + end subroutine calculate_volume + + +end module density diff --git a/src/data_input.f90 b/src/input_data/data_input.f90 similarity index 76% rename from src/data_input.f90 rename to src/input_data/data_input.f90 index 4179d92..a0b658a 100644 --- a/src/data_input.f90 +++ b/src/input_data/data_input.f90 @@ -42,18 +42,22 @@ subroutine read_components(file, def_components, scn, scn_plus) close(funit) end subroutine read_components - subroutine read_composition(file, def_comp_z, scn_z, plus_z) + subroutine read_composition(file, def_comp_z, scn_z, plus_z, plus6_z_exp, plus7_z_exp, plus12_z_exp, plus30_z_exp ) !! Reads the molar compositions of each component from input file character(len=*), intent(in) :: file !! file name real(pr), allocatable, intent(out) :: def_comp_z(:) !! set of corresponding mole fractions of defined components real(pr), allocatable, intent(out) :: scn_z(:) !! set of corresponding mole fractions of scn cuts real(pr), intent(out):: plus_z !! composition of residual fraction from input file + real(pr), intent(out):: plus6_z_exp !! composition of residual fraction from input file + real(pr), intent(out):: plus7_z_exp !! composition of residual fraction from input file + real(pr), intent(out):: plus12_z_exp !! composition of residual fraction from input file + real(pr), intent(out):: plus30_z_exp !! composition of residual fraction from input file integer :: def_comp_nc !! number of defined components being considered in the oil integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped integer :: funit - namelist /nml_composition/ def_comp_z, scn_z, plus_z + namelist /nml_composition/ def_comp_z, scn_z, plus_z, plus6_z_exp, plus7_z_exp, plus12_z_exp, plus30_z_exp call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(def_comp_z(def_comp_nc)) @@ -64,18 +68,22 @@ subroutine read_composition(file, def_comp_z, scn_z, plus_z) close(funit) end subroutine read_composition - subroutine read_molecular_weight(file, def_comp_mw, scn_mw, plus_mw) + subroutine read_molecular_weight(file, def_comp_mw, scn_mw, plus_mw, plus6_mw_exp, plus7_mw_exp, plus12_mw_exp, plus30_mw_exp) !! Reads the molecular weights of each component from the input file character(len=*), intent(in) :: file !! file name real(pr), allocatable, intent(out) :: def_comp_mw(:) !! set of corresponding molecular weights of defined components real(pr), allocatable, intent(out) :: scn_mw(:) !! set of corresponding molecular weights of scn cuts real(pr), intent(out) :: plus_mw !! molecular weight of residual fraction + real(pr), intent(out) :: plus6_mw_exp !! molecular weight of residual fraction + real(pr), intent(out) :: plus7_mw_exp !! molecular weight of residual fraction + real(pr), intent(out) :: plus12_mw_exp !! molecular weight of residual fraction + real(pr), intent(out) :: plus30_mw_exp !! molecular weight of residual fraction integer :: def_comp_nc !! number of defined components being considered in the oil integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped integer :: funit - namelist /nml_molecular_weight/ def_comp_mw, scn_mw, plus_mw + namelist /nml_molecular_weight/ def_comp_mw, scn_mw, plus_mw, plus6_mw_exp, plus7_mw_exp, plus12_mw_exp, plus30_mw_exp call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(def_comp_mw(def_comp_nc)) @@ -104,9 +112,17 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts real(pr) :: plus_z !! composition of residual fraction + real(pr) :: plus6_z_exp !! composition of residual fraction from input file + real(pr) :: plus7_z_exp !! composition of residual fraction from input file + real(pr) :: plus12_z_exp !! composition of residual fraction from input file + real(pr) :: plus30_z_exp !! composition of residual fraction from input file real(pr), allocatable :: def_comp_mw(:) !! set of corresponding molecular weights of defined components real(pr), allocatable :: scn_mw(:) !! set of corresponding molecular weights of scn cuts real(pr) :: plus_mw !! molecular weight of residual fraction + real(pr) :: plus6_mw_exp !! molecular weight of residual fraction + real(pr) :: plus7_mw_exp !! molecular weight of residual fraction + real(pr) :: plus12_mw_exp !! molecular weight of residual fraction + real(pr) :: plus30_mw_exp !! molecular weight of residual fraction real(pr), allocatable :: product_z_mw_i(:) real(pr), allocatable, intent(out) :: def_comp_w(:) real(pr), allocatable, intent(out) :: scn_w(:) @@ -124,8 +140,12 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& allocate(def_comp_w(def_comp_nc)) allocate(scn_w(scn_nc)) - call read_composition(file, def_comp_z, scn_z, plus_z) - call read_molecular_weight(file, def_comp_mw, scn_mw, plus_mw) + call read_composition(file, def_comp_z, scn_z, plus_z, plus6_z_exp, & + plus7_z_exp, plus12_z_exp, plus30_z_exp & + ) + call read_molecular_weight(file, def_comp_mw, scn_mw, plus_mw, & + plus6_mw_exp, plus7_mw_exp, plus12_mw_exp, plus30_mw_exp & + ) product_z_mw_def_comp = (def_comp_z)*(def_comp_mw) product_z_mw_scn = (scn_z)*(scn_mw) @@ -138,21 +158,27 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& plus_w = w(def_comp_nc+scn_nc+1) end subroutine mass_fractions - subroutine read_density(file,scn_density,plus_density) + subroutine read_density(file,scn_density,plus_density, plus6_density_exp, & + plus7_density_exp, plus12_density_exp, plus30_density_exp & + ) !! Reads the density of each component from the input file and calculated molar volume since C6 fraction. character(len=*), intent(in) :: file !! file name real(pr), allocatable, intent(out) :: scn_density(:) !! set of corresponding densities of scn cuts real(pr), intent(out):: plus_density !! density of residual fraction from input file + real(pr), intent(out):: plus6_density_exp + real(pr), intent(out):: plus7_density_exp + real(pr), intent(out):: plus12_density_exp + real(pr), intent(out):: plus30_density_exp integer :: def_comp_nc !! number of defined components being considered in the oil integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped integer :: funit - namelist /nml_density/ scn_density, plus_density - + namelist /nml_density/ scn_density, plus_density, plus6_density_exp, plus7_density_exp, plus12_density_exp, plus30_density_exp + call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(scn_density(scn_nc)) - + open(newunit=funit, file=file) read(funit, nml=nml_density) close(funit) @@ -162,22 +188,31 @@ end subroutine read_density type(FluidData) function data_from_file(file) result(data) !! This funtion allows to obtain experimental data from data imput character(len=*), intent(in) :: file !! file name - + data%filename = file - - call read_setup(file, data%def_comp_nc, data%scn_nc, data%scn_nc_ps, data%numbers_ps) - call read_components(file, data%def_components, data%scn, & - data%scn_plus) - call read_composition(file, data%def_comp_z, data%scn_z, data%plus_z) + + call read_setup(file, data%def_comp_nc, data%scn_nc, data%scn_nc_ps, & + data%numbers_ps & + ) + call read_components(file, data%def_components, data%scn, data%scn_plus) + call read_composition(file, data%def_comp_z, data%scn_z, data%plus_z, & + data%plus6_z_exp, data%plus7_z_exp, data%plus12_z_exp, data%plus30_z_exp & + ) call read_molecular_weight(file, data%def_comp_mw, data%scn_mw, & - data%plus_mw) + data%plus_mw, data%plus6_mw_exp, data%plus7_mw_exp, data%plus12_mw_exp, & + data%plus30_mw_exp & + ) call mass_fractions(file, data%w, data%product_z_mw_def_comp, & data%product_z_mw_scn, data%sum_z_mw_i, data%product_z_mw_plus, & - data%def_comp_w, data%scn_w, data%plus_w) - call read_density(file, data%scn_density, data%plus_density) - + data%def_comp_w, data%scn_w, data%plus_w & + ) + call read_density(file, data%scn_density, data%plus_density, & + data%plus6_density_exp, data%plus7_density_exp, data%plus12_density_exp, & + data%plus30_density_exp & + ) + end function data_from_file - + end module data_from_input diff --git a/src/lumping/lumping_method.f90 b/src/lumping/lumping_method.f90 new file mode 100644 index 0000000..e082987 --- /dev/null +++ b/src/lumping/lumping_method.f90 @@ -0,0 +1,184 @@ +module lumping + use dtypes + +contains + + subroutine lump(fluid, characterization) + !! This sobroutine ... + implicit none + type(FluidData), intent(in) :: fluid + type(FluidDataOut), intent(inout) :: characterization + integer :: last_C ! CN of last single cut: typically 19 + integer :: i_last ! Number of elements in the distribution of CN+ + integer :: last ! Total elements starting after defined components. + integer :: scn_nc_input !! inicial number of single cuts being considered in the oil from data input + integer :: scn_nc_new !! new number of single cuts being considered in the oil defined as from scn_nc_ps variable. + real(pr) :: plus_w_new + real(pr), allocatable :: z_m_plus_i(:) + real(pr), allocatable :: plus_z_i(:) + real(pr), allocatable :: var_aux_1(:) + real(pr), allocatable :: var_aux_2(:) + integer :: i, prev_i + integer, dimension(15) :: j_ps + integer :: i_ps ! auxiliary variables + real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps ! auxiliary variables + real(pr), allocatable :: plus_z_ps(:) + real(pr), allocatable :: plus_mw_ps(:) + integer :: numbers_ps + real(pr), allocatable :: density_ps(:) + real(pr), allocatable :: w_ps(:) + real(pr) :: sum_zm_last_ps + real(pr), allocatable :: moles(:) + real(pr), allocatable :: lumped_z(:) + real(pr), allocatable :: lumped_mw(:) + real(pr), allocatable :: lumped_densities(:) + + associate (& + def_nc => fluid%def_comp_nc , & + scn_mw => characterization%scn_mw, & + plus_zm => characterization%plus_zm, & + plus_density => characterization%plus_density, & + scn_zm => characterization%scn_zm, & + plus_z => characterization%plus_z, & + scn_z => characterization%scn_z, & + plus_w => fluid%plus_w, & + w => fluid%w, & + plus_mw => characterization%plus_mw, & + carbon_number_plus => characterization%carbon_number_plus, & + plus6_density => characterization%plus6_density, & + C => characterization%C, & + scn_nc_new => characterization%scn_nc_new & + ) + + last_C = fluid%scn(fluid%scn_nc) !19 + i_last = characterization%nc_plus !124 + last = fluid%scn_nc + i_last !138 + + allocate (z_m_plus_i(0)) + allocate (plus_z_i(0)) + allocate (var_aux_1(0)) + allocate (var_aux_2(0)) + allocate (plus_z_ps(0)) + allocate (plus_mw_ps(0)) + allocate(density_ps(0)) + allocate (w_ps(0)) + allocate(moles(0)) + allocate(lumped_z(0)) + allocate(lumped_mw(0)) + allocate(lumped_densities(0)) + + + scn_nc_input = fluid%scn_nc + scn_nc_new = fluid%scn_nc_ps - 6 + ! scn_nc_ps : CN from which all SCN fractions will be lumped + ! into the specified number of pseudos + + characterization%plus_w = plus_w + z_m_plus_i = characterization%product_z_mw_plus_i + plus_z_i = characterization%plus_z_i + + if (scn_nc_new < scn_nc_input) then + ! Plus Fraction needs to be extended to include lower CN's + plus_zm = plus_zm + sum(scn_zm(scn_nc_new + 1 : scn_nc_input)) + plus_z = plus_z + sum(scn_z(scn_nc_new + 1 : scn_nc_input)) !extended zp + plus_w_new = plus_w + sum(w(fluid%def_comp_nc+scn_nc_new + 1 & + : fluid%def_comp_nc+scn_nc_input)) + !! use new variable for plus_w because type fluidata can't be modified. + plus_mw = plus_zm / plus_z + z_m_plus_i = scn_zm(scn_nc_new + 1 : scn_nc_input ) + z_m_plus_i = [z_m_plus_i, characterization%product_z_mw_plus_i] + plus_z_i = scn_z(scn_nc_new + 1 : scn_nc_input ) + plus_z_i = [plus_z_i, characterization%plus_z_i] + i_last = i_last + scn_nc_input - scn_nc_new + last_C = last_C - scn_nc_input + scn_nc_new + characterization%plus_w = plus_w_new + characterization%product_z_mw_plus_i = z_m_plus_i + characterization%plus_z_i = plus_z_i + end if + + !do i= 1, scn_nc_new + ! print*, scn_z(i), scn_mw(i) + !end do + + numbers_ps = fluid%numbers_ps + j_ps = 0 + ! Lumping into Nps pseudos + rec_zm = plus_zm / numbers_ps + ! Recommended value for the product z*M (proportional to weight) + ! for each pseudo according to pedersen + remain_plus_zm = plus_zm + i_ps = 1._pr + sum_z = 0.0_pr + sum_zm = 0.0_pr + sum_volume_ps = 0.0_pr + + i = 0.0_pr + + do while (i_ps < numbers_ps) + i = i + 1 + j_ps(i_ps) = j_ps(i_ps) + 1 + sum_z = sum_z + plus_z_i(i) + sum_zm = sum_zm + z_m_plus_i(i) + sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / plus6_density(scn_nc_new+i) + if (z_m_plus_i(i+1) > 2*(rec_zm - sum_zm)) then + ! when adding one more would go too far + plus_z_ps = [plus_z_ps, sum_z] + plus_mw_ps = [plus_mw_ps, sum_zm/sum_z] + w_ps = [w_ps, (characterization%plus_w*(sum_zm/plus_zm))] + remain_plus_zm = remain_plus_zm - sum_zm + if (remain_plus_zm < rec_zm) numbers_ps = i_ps + 1 + carbon_number_plus(scn_nc_new + i_ps) = 6+((plus_mw_ps(i_ps)-84)/C) + density_ps = [density_ps, (sum_zm/sum_volume_ps)] + i_ps = i_ps + 1 + sum_z = 0._pr + sum_zm = 0._pr + sum_volume_ps = 0._pr + end if + end do + + plus_z_ps = [plus_z_ps, sum(plus_z_i(i+1:i_last))] + ! at this point, Nps is the order for the last NON-Asphaltenic pseudo comp. + ! (e.g. 4 if 5 is for Asphaltenes) + plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i+1:i_last))/plus_z_ps(numbers_ps)] + w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps))* & + (plus_mw_ps(numbers_ps))/(plus_zm))] + carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & + ((plus_mw_ps(numbers_ps) - 84) / C) + sum_zm_last_ps = (sum(z_m_plus_i(i+1:i_last))) + density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / (plus6_density(scn_nc_new+i+1:last)))] + j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) !! revisar manana + + !do i = 1, numbers_ps + ! print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) + !end do + + plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & + plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) + + moles = [(fluid%def_comp_w)/(fluid%def_comp_mw)] + moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] + moles = [moles, w_ps/plus_mw_ps ] + + characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize + characterization%last_C = last_C + characterization%i_last = i_last + characterization%last = last + + ! In this point we create a new array to content final Mis and densities + ! whit pseudos. + + lumped_z = [scn_z(1:scn_nc_new), plus_z_ps] + lumped_mw = [scn_mw(1:scn_nc_new), plus_mw_ps] + lumped_densities = [plus6_density(1:scn_nc_new), density_ps ] + + + characterization%lumped_z = lumped_z + characterization%lumped_mw = lumped_mw + characterization%lumped_densities = lumped_densities + + end associate + + end subroutine lump + +end module lumping + diff --git a/src/module.f90 b/src/module.f90 deleted file mode 100644 index 19e8c82..0000000 --- a/src/module.f90 +++ /dev/null @@ -1,836 +0,0 @@ - -module data - - use constants, only: pr - - !use types, only: FluidData - - implicit none - integer :: nv, nv1, nNcut, nNdef - real(pr) :: Zplus - - -contains - - subroutine data_in (Nfluid,Oil,Ncut,Ndef,DefComp,zdef,rMdef,rCN,zcomp,rMW,& - Den,Plus, rMWplus, DenPlus,Nps,wat,watPlus,zM,zMp,Z6p,sumV,w, rn ) - - implicit none - - integer, parameter :: maxD=15, imax=48 - integer :: i,j!, nv1 - integer :: file_unit_1, file_unit_2, file_unit_3, Nfluid, Ncut, Ndef, Nps - integer, dimension(330) :: rCN - real(pr), dimension(maxD) :: zdef, rMdef - real(pr), dimension(imax) :: zcomp, rMW, Den, wat, zM, zMdef, w, rn,MWt - real(pr) :: rMWplus, DenPlus , watPlus, zMp, Z6p, sumV, zMtotal - character(len=30) :: infile, outfile - character(len=7) :: Oil, Plus - character(len=4) :: DefComp(maxD) - - write(*,*) 'ARE MOLECULAR WEIGTHS FROM DATAIN EXPERIMENTAL?' - write(*,*) '1: YES 2: NO' - read(*,*) nv - - write(*,*) 'ENTER THE VERSION NUMBER' - write(*,*) '1: C.M 2: CMplus' - read(*,*) nv1 - - !nv2 = nv1 - - ! open input file - write(*,*) 'ENTER INFILE' - read(*,*) infile - write(*, *) "READING INFILE: " // infile - open(newunit=file_unit_1, file =infile) - write(*,*) 'ENTER OUTFILE' - read(*,*) outfile - open(newunit=file_unit_2, file = outfile) - open(newunit=file_unit_3, file = 'TfPsatIn.txt') - ! read input file:read number fluid, single carbon number for each fluid. - read(file_unit_1,*) Nfluid - do j = 1, Nfluid - read (file_unit_1,*) Oil - read (file_unit_1,*) Ncut, Ndef - end do - - nNcut = Ncut - nNdef = Ndef - - - ! read defined compounds with molar fraction and M by fluid j and compound i. - do i = 1, Ndef - read(file_unit_1,*) DefComp(i), zdef(i), rMdef(i) - end do - ! read carbon number and molar fraction of fluid j and comp. i - do i = 1, Ncut - read(file_unit_1,*) rCN(i), zcomp(i), rMW(i),Den(i) - end do - - read(file_unit_1,*) Plus, zPlus, rMWplus, DenPlus ! plus fraction properties - read (file_unit_1,*) Nps !number of pseudos for the C20+ fraction - read (file_unit_1,*) wat(1:Ndef+Ncut),watPlus ! atmospheric oil weight fractions - - !zpl = zPlus ! added by oscar - ! There get Zi*Mi and volume from cut to plus fraction - - zMdef(1:Ndef) = zdef(1:Ndef)*rMdef(1:Ndef) - zM(1:Ncut) = zcomp(1:Ncut) * rMW(1:Ncut) - zMp = zPlus * rMWplus - zMtotal = sum(zMdef(1:Ndef)) + sum(zM(1:Ncut)) + zMp - w(1:Ndef) = zMdef(1:Ndef) / zMtotal - w(Ndef+1:Ndef+Ncut) = zM(1:Ncut) / zMtotal - w(Ndef+Ncut+1) = zMp / zMtotal - rn(1:Ndef) = zdef(1:Ndef)/zMtotal - Z6p = sum(zcomp(1:Ncut)) + zPlus - sumV = sum(zM(1:Ncut)/Den(1:Ncut)) + zMp/DenPlus - - !agregado para calcular el peso molecular global - !MWt(1:Ndef) = rMdef(1:Ndef) - !MWt(Ndef+1:Ndef+Ncut) = rMW(1:Ncut) - - !do i= 1, Ndef+Ncut - ! print*, MWt(i) - !end do - - - - - !print*, Zplus - end subroutine data_in - -end module data - -module routines_pre - - use constants, only: pr - !use dtypes, only: FluidDataOut - use data - - - - contains - - subroutine Linear_Regression(x,y,a,b,r2) - !! This subroutine computes the regression line for a data set of x, y variables. - implicit none - - integer :: i - integer :: n !! total number of data set. - real(pr), allocatable, intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable - real(pr), allocatable, intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable - real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables - real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line - real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line - real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient - - n = size(x) - !Calculation of a y b --> y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 - end do - - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 - end do - - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) - - end subroutine Linear_Regression - - - subroutine Best_Linear_Regression(scn_nc,scn,scn_z,plus_z,a,b,r2,ninit) - !! This subroutine calculates the best regression line for an oil. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer :: i, j, k , kold, nbest, xaux,cbmax ! internal variables - real, allocatable, intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: ninit !! minimum carbon number obtained from the best linear regression - real(pr), allocatable :: xBR(:), yBR(:) - real(pr) :: r2old,r2best, aold, bold, abest, bbest, zsum,zaux - - allocate(scn(scn_nc)) - allocate(scn_z(scn_nc)) - allocate(xBR(scn_nc)) - allocate(yBR(scn_nc)) - - k=5 - r2=0.0001d0 - r2old=0.00001d0 - r2best=0.0001d0 - - do while (r2.gt.r2old.or.r2old.lt.0.9) - kold=k - r2old=r2 - aold=a - bold=b - - if (r2.gt.r2best) then - r2best=r2 - abest=a - bbest=b - nbest=scn(scn_nc-k+2) - end if - - if (k.gt.scn_nc) then - if (r2.gt.r2best)then - ninit=scn(scn_nc-k+2) - go to 22 - else - r2=r2best - a=abest - b=bbest - ninit=nbest - go to 22 - end if - end if - - j=1 - xBR=0. - yBR=0. - - do i=scn_nc-k+1, scn_nc - xBR(j)=scn(i) - yBR(j)=scn_z(i) - j=j+1 - end do - - call LinearRegression(xBR(:k), yBR(:k), a, b, r2) - k=k+1 - end do - - r2=r2old - a=aold - b=bold - ninit=scn(scn_nc-kold+2) - - 22 continue - - zsum = 0d0 - xaux=scn(scn_nc) - - do while (zsum.lt.plus_z.and.xaux<300) - xaux=xaux+1 - zaux= exp(a*xaux+b) - zsum=zsum+zaux - end do - - cbmax=xaux - - print*, a, b, cbmax, r2, Ninit - - end subroutine Best_Linear_Regression - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - subroutine get_C(start,C,Ncut,rCN,zM,zMp,Z6p,a,b,rMp,maxC,z,zpi,zMpi,Ndef,zp,w,rn,rMWplus,zcomp,zdef) - - implicit none - - integer, parameter :: imax=48, maxD=15 - logical :: start - real(pr) :: C, Z6p, zMp, a, b, rMp, rMWplus - real(pr) :: aux, dold, dif, Cold,zp - real(pr), dimension(imax) :: zM ! inputs coming from main program - integer, dimension(imax) :: rCN ! inputs coming from main program - real(pr), dimension(imax) :: z, rM,w,rn, zcomp - integer :: Ncut, maxC, Ndef - real(pr), dimension(300):: zpi, zMpi - real(pr), dimension(maxD) :: zdef - - !write(*,*) 'ENTER THE VERSION NUMBER' - !write(*,*) '1: CMOD 2: CMM' - !read(*,*) nv - - - if(nv1==1)then - - ! Find the C value for molecular weights line - C = 13.0 ! initial guess - Start = .true. - rMp = rMWplus - - call difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dold,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - Cold = C - C = min(12.8,C-0.07) - Start = .false. - - - - call difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dif,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - - do while (abs(dif) > 0.1) - aux = C - C = C - dif*(C-Cold)/(dif-dold) - Cold = aux - dold = dif - call difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dif,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - end do - - !print*, C,maxC, rMp - - - else - - ! Find the C value for molecular weights line - C = 14.0 ! initial guess and maximum C allowed - Start = .true. - rMp = rMWplus - - call difMpfromCfull (Start,C,Ndef,Ncut,rCN,dold,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef) - Cold = C - C = 13.5 !min(12.8,C-0.07) - Start = .false. - - call difMpfromCfull (Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn, zcomp,zdef) - - do while (abs(dif) > 0.1) - aux = C - C = C - dif*(C-Cold)/(dif-dold) - Cold = aux - dold = dif - call difMpfromCfull (Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef) - end do - - !print*, C, maxC, rMp, zp - - end if - - end subroutine get_C - - subroutine difMpfromC(start,C,Ncut,rCN,zM,zMp,Z6p,dif,a,b,rMp,maxC,z,zpi,zMpi,zcomp,zdef,zp) - - use data - - implicit none - - integer, parameter :: imax=48, maxD=15 - logical :: start - real(pr) :: C, Z6p, zMp, dif, a, b, rMp, zp,aBE,bBE,r2, alim, blim, half - real(pr) :: a60,b60, sumz, rMpcalc - real(pr), dimension(imax) :: zM ! inputs coming from main program - integer, dimension(imax) :: rCN,x ! inputs coming from main program - real(pr), dimension(imax) :: ylog, rM, z, zcomp - integer :: i, Ncut, maxC, i0, CmaxL, C60max - real(pr), dimension(300):: zpi, zMpi - real(pr), dimension(maxD) :: zdef - - ! 1 i0 = rCN(Ncut) - - i = 300 - sumz = -1 - do while (i==300.and.sumz12 y C<12. - a = alim - b = blim - else - !call LineC60max (Ncut,Zp,rCN,aBE,bBE,half,a60,b60,C60max) - if(aBE > a60)then - a = a60 - b = b60 - else - a = aBE - b = bBE - end if - end if - - sumz = 0.0d0 - i = 0 - do while (sumz a60)then - a = a60 - b = b60 - else - a = aBE - b = bBE - end if - end if - - - sumz = 0.0d0 - i = 0 - do while (sumzalim) eliminado por oscar 20/03/2024 - ! aold = a - ! a = max(1.0001*a, alim) - ! b = b - (a-aold)*19.51 - ! sumz = 0.0d0 - ! i = 0 - ! do while (sumz 0.00001) - aux = rMp - rMp = rMp - dif*(rMp-rMpold)/(dif-dold) - rMpold = aux - dold = dif - call difnewMp(Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef,rMW,rMdef, Mglobal) - end do - - - - - !print*, C, maxC, rMp - !print*, zp , z6p - !print*, z, zp - end if - - - end subroutine GetNewMpfromC - - - subroutine difnewMp(Start,C,Ndef,Ncut,rCN,dif,a,b,rM,rMp,maxC,z,zp,zpi,zMpi,w,rn,zcomp,zdef,rMW,rMdef, Mglobal) - !! this subroutines for a C value returns the correponding Mp value - !use MassFracAndMoles, only: w,rn - use data - implicit none - - integer, parameter :: imax=48, maxD=15 - logical :: start - real(pr) :: C, dif, a, b, rMp, zp ,aBE ,bBE ,r2 , alim, blim, half - real(pr) :: a60,b60, sumz, rMpcalc, rntot, aold, Mglobal - integer, dimension(imax) :: rCN,x ! inputs coming from main program - real(pr), dimension(imax) :: ylog, rM, z, w, rn, zcomp, ztotal,rMW, MWt - integer :: i,Ndef, Ncut, maxC, i0, CmaxL, C60max - real(pr), dimension(300):: zpi, zMpi,ZM_1 - real(pr), dimension(maxD) :: zdef,rMdef - - 1 i0 = rCN(Ncut) - - if(nv == 1)then !introduce 'if' to choise where experimental data from - x = rCN - rn(Ndef+1:Ndef+Ncut) = w(Ndef+1:Ndef+Ncut)/rMW(1:Ncut) - rn(Ndef+Ncut+1) = w(Ndef+Ncut+1)/rMp - rntot = sum(rn(1:Ndef+Ncut+1)) - ztotal(1:Ndef+Ncut) = rn(1:Ndef+Ncut) / rntot - z(1:Ncut) = rn(Ndef+1:Ndef+Ncut) / rntot - Zp = rn(Ndef+Ncut+1) / rntot - !ylog = log(zcomp) - ylog = log(z) - - !calculo m global - MWt(1:Ndef) = rMdef(1:Ndef) - MWt(Ndef+1:Ndef+Ncut) = rMW(1:Ncut) - ZM_1(1:Ndef+Ncut) = ztotal(1:Ndef+Ncut)*MWt(1:Ndef+Ncut) - - !do i= 1, Ndef+Ncut - ! print*, MWt(i) - !end do - - end if - - if(nv == 2)then - rM(1:Ncut) = 84 + C*(rCN(1:Ncut)-6) - rn(Ndef+1:Ndef+Ncut) = w(Ndef+1:Ndef+Ncut)/rM(1:Ncut) - rn(Ndef+Ncut+1) = w(Ndef+Ncut+1)/rMp - rntot = sum(rn(1:Ndef+Ncut+1)) - ztotal(1:Ndef+Ncut) = rn(1:Ndef+Ncut) / rntot - z(1:Ncut) = rn(Ndef+1:Ndef+Ncut) / rntot - Zp = rn(Ndef+Ncut+1) / rntot - x = rCN - ylog = log(z) - !calculo m global - MWt(1:Ndef) = rMdef(1:Ndef) - MWt(Ndef+1:Ndef+Ncut) = rM(1:Ncut) - ZM_1(1:Ndef+Ncut) = ztotal(1:Ndef+Ncut)*MWt(1:Ndef+Ncut) - - !do i= 1, Ndef+Ncut - ! print*, MWt(i) - !end do - - end if - !print*, zp - !print*, ztotal - - call BestLinearRegression(Ncut,x,ylog,Zp,aBE,bBE,r2) - call LimitLine(Ncut,rCN,Zp,aBE,bBE,half,alim,blim,CmaxL) - call LineC60max (Ncut,Zp,rCN,aBE,bBE,half,a60,b60,C60max) - - if(aBE < alim)then !, para C se excede el valor de 14 en la mayoria de los casos. - a = alim - b = blim - else - !call LineC60max (Ncut,Zp,rCN,aBE,bBE,half,a60,b60,C60max) - if(aBE > a60)then - a = a60 - b = b60 - else - a = aBE - b = bBE - end if - end if - - - sumz = 0.0d0 - i = 0 - do while (sumzalim) - aold = a - a = max(1.0001*a, alim) - b = b - (a-aold)*19.51 - sumz = 0.0d0 - i = 0 - do while (sumz fluid%scn_nc, def_comp_nc => fluid%def_comp_nc, & + def_comp_w => fluid%def_comp_w, def_comp_mw => fluid%def_comp_mw, & + scn_w => fluid%scn_w, plus_w => fluid%plus_w, scn => fluid%scn, & + plus_mw => characterization%plus_mw, scn_z => characterization%scn_z, & + plus_z => characterization%plus_z, log_scn_z => characterization%log_scn_z, & + C => characterization%C, scn_mw => characterization%scn_mw, & + scn_zm => characterization%scn_zm, & + plus_zm => characterization%plus_zm & + ) + + allocate(def_comp_moles(def_comp_nc)) + allocate(scn_moles(scn_nc)) + + def_comp_moles = (def_comp_w)/(def_comp_mw) + + select case (mw_source) + + case("experimental") + scn_mw = fluid%scn_mw + scn_moles = (scn_w)/(scn_mw) + plus_moles = (plus_w)/(plus_mw) + total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) + scn_z = scn_moles / total_moles + plus_z = plus_moles / total_moles + + case("calculated") + + if (method=="global_mw")then + scn_mw = 84 + C*(scn-6) + scn_z = fluid%product_z_mw_scn / scn_mw + sum_def_comp_z_plus = sum(fluid%scn_z) + (fluid%plus_z) + plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C + plus_mw = fluid%product_z_mw_plus / plus_z + end if + + if (method=="plus_mw")then + scn_mw = 84 + C*(scn-6) + scn_moles = (scn_w)/(scn_mw) + plus_moles = (plus_w)/(plus_mw) + total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) + scn_z = scn_moles / total_moles + plus_z = plus_moles / total_moles + + endif + end select + scn_zm = scn_z*scn_mw + plus_zm = plus_z*plus_mw + log_scn_z = log(scn_z) + end associate + + end subroutine select_method + + subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & + characterization) + + !! this subroutine ... + implicit none + + type(FluidData), intent(inout) :: fluid + type(FluidDataOut), intent(inout) :: characterization + logical :: start + character(len=*), intent(in), optional :: method + character(len=*), intent(in) :: mw_source + real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction + real(pr), intent(out) :: difference + real(pr) :: half + real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux + real(pr):: sum_z + real(pr) :: denom + !! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] + real(pr) :: sum_def_comp_z_plus + !! compositions sum from (define component +1) to (plus component) + integer :: j, k , k_old, n_best, x_aux, i_0! internal variables + integer :: i + integer, dimension(300) :: carbon_number_plus + real(pr), dimension(300) :: plus_z_i + real(pr), dimension(300) :: product_z_mw_plus_i + real(pr), dimension(300) :: scn_i + + associate (& + scn_nc => fluid%scn_nc, scn => fluid%scn,& + plus_z => characterization%plus_z, plus_mw => characterization%plus_mw, & + log_scn_z => characterization%log_scn_z, C => characterization%C, & + a_blr => characterization%a_blr, b_blr => characterization%b_blr, & + a_60 => characterization%a_60, b_60 => characterization%b_60, & + a_lim => characterization%a_lim, b_lim => characterization%b_lim, & + r2 => characterization%r2, n_init => characterization%n_init, & + c_max_blr => characterization%c_max_lim, a => characterization%a, & + c_max_lim => characterization%c_max_lim, b => characterization%b, & + c_max_60 => characterization%c_max_lim, & + c_max => characterization%c_max & + ) + + 1 i_0 = scn(scn_nc) + + call select_method(fluid, mw_source, method, characterization) ! select case + call Best_Linear_Regression(scn_nc, scn, log_scn_z, & + plus_z, a_blr, b_blr, r2, n_init, c_max_blr) + call LimitLine(scn_nc, scn, plus_z, a_blr, b_blr,a_lim,b_lim,c_max_lim,half) + call Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. + + if(a_blr < a_lim)then + !Added by oscar 05/12/2023.se elimina por que se agregan las restricciones + !para characterization%C>12 y characterization%C<12. + a = a_lim + b = b_lim + else ! this line was removed call Line_C60_max + if(a_blr > a_60)then + a = a_60 + b = b_60 + else + a = a_blr + b = b_blr + end if + end if + + sum_z = 0.0_pr + i = 0 + do while (sum_z < plus_z .and. i< 300) + i = i+1 + plus_z_i(i) = exp(a*(i+i_0)+b) + sum_z = sum_z + plus_z_i(i) + product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) + carbon_number_plus(i) = i+i_0 + end do + + if (i == 300.and.sum_z < plus_z)then + if(start) C = C - 0.07 + if(.not.start) C = C - 0.01 + go to 1 + end if + + c_max = i+i_0 + plus_z_i(i) = plus_z_i(i) - (sum_z-plus_z) ! Adjustment to Zp (z20+) + product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) + + if (mw_source == "experimental" .and. C>12 .and. C<14 ) then + ! mayor igual a 12 o menor igual 14. + denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) + ! denominator of equation to obtain C value directly + C = (plus_z*(plus_mw-84))/denom + endif + + plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z + difference = plus_mw_cal-plus_mw + scn_i = [scn, carbon_number_plus(1:i)] + end associate + + characterization%plus_z_i = plus_z_i(1:i) + characterization%product_z_mw_plus_i = product_z_mw_plus_i(1:i) + characterization%carbon_number_plus = carbon_number_plus(1:i) + characterization%nc_plus = i + characterization%scn_i = scn_i + + + + end subroutine difference_mw_plus + + subroutine get_c_or_m_plus(fluid, mw_source, method, fix_C, characterization) + !! This subroutine... + implicit none + type(FluidData) :: fluid + type(FluidDataOut) :: characterization + logical :: start + character(len=*), intent(in) :: mw_source + character(len=*), intent(in), optional :: method + logical, intent(in) :: fix_C !! If C is < 12 or > 14 then fix to the + !! closest value. + !integer :: c_max !! output CN at which characterization%plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) + real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cutsn + real(pr) :: difference !! + real(pr) :: difference_old, plus_mw_old + real(pr) :: C_old + real(pr) :: aux + + associate (C => characterization%C, plus_mw => characterization%plus_mw) + + select case (mw_source) + case("experimental") + C = 13.5 + Start = .true. + plus_mw = fluid%plus_mw + call difference_mw_plus(fluid, mw_source, method, start, difference, & + characterization) + + case("calculated") + + if(method == 'global_mw') C=13 + if(method == 'plus_mw') C=14 + + Start = .true. + plus_mw = fluid%plus_mw + call difference_mw_plus(fluid, mw_source, method, start, & + difference_old, characterization) + C_old = C + if(method == 'global_mw') C = min(12.8, C-0.07) + if(method == 'plus_mw') C = 13.5 + Start = .false. + call difference_mw_plus(fluid, mw_source, method, start, difference, & + characterization) + + do while (abs(difference) > 0.1) + aux = C + C = C - difference*(C-C_old)/(difference-difference_old) + C_old = aux + difference_old = difference + call difference_mw_plus(fluid, mw_source, method, start, & + difference, characterization) + end do + end select + + if(fix_C .and. C>14 .or. C<12)then + + + if(C>14) C=14 + if(C<12) C=12 + + Start = .true. + plus_mw = fluid%plus_mw + call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & + difference_old, characterization) + plus_mw_old = plus_mw + plus_mw = 0.9*plus_mw + Start = .false. + call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & + difference, characterization) + + do while (abs(difference)>0.00001) + aux = plus_mw + plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) + plus_mw_old = aux + difference_old = difference + call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & + difference, characterization) + end do + end if + end associate + end subroutine get_c_or_m_plus + +end module molecular_weigth + diff --git a/src/prueba.f90 b/src/prueba.f90 deleted file mode 100644 index 0a787d7..0000000 --- a/src/prueba.f90 +++ /dev/null @@ -1,466 +0,0 @@ -module pruebas - use constants - use data_from_input, only: data_from_file, FluidData - - - contains - - subroutine Linear_Regression(x,y,a,b,r2) - !! This subroutine computes the regression line for a data set of x, y variables. - implicit none - - integer :: i - integer :: n !! total number of data set. - real(pr), intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable - real(pr), intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable - real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables - real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line - real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line - real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient - - n = size(x) - !Calculation of a y b --> y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 - end do - - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 - end do - - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) - - end subroutine Linear_Regression - - subroutine Best_Linear_Regression(scn_nc,scn,scn_z,plus_z,a,b,r2,n_init,c_max_blr) - !! This subroutine calculates the best regression line for an oil. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer, intent(out) :: c_max_blr - !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable, intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression - integer :: i, j, k , k_old, n_best, x_aux ! internal variables - real(pr), dimension(scn_nc) :: x_blr, y_blr - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - - k=5 - r2=0.0001_pr - r2_old=0.00001_pr - r2_best=0.0001_pr - - do while (r2.gt.r2_old.or.r2_old.lt.0.9) - k_old=k - r2_old=r2 - a_old=a - b_old=b - - if (r2.gt.r2_best) then - r2_best=r2 - a_best=a - b_best=b - n_best=scn(scn_nc-k+2) - end if - - if (k.gt.scn_nc) then - if (r2.gt.r2_best)then - n_init=scn(scn_nc-k+2) - go to 22 - else - r2=r2_best - a=a_best - b=b_best - n_init=n_best - go to 22 - end if - end if - - j=1 - x_blr=0. - y_blr=0. - - do i=scn_nc-k+1, scn_nc - x_blr(j)=scn(i) - y_blr(j)=scn_z(i) - j=j+1 - end do - - call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) - k=k+1 - end do - - r2=r2_old - a=a_old - b=b_old - n_init=scn(scn_nc-k_old+2) - - 22 continue - - z_sum = 0d0 - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a*x_aux+b) - z_sum=z_sum+z_aux - end do - - c_max_blr=x_aux - - !print*, a, b, c_max_blr, r2, n_init - - end subroutine Best_Linear_Regression - - subroutine LimitLine(scn_nc,scn,plus_z,a_blr,b_blr,a_lim,b_lim,c_max_lim,half) - !! This subroutine obtains the limit line constants for an oil. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer, intent(out) :: c_max_lim - !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr) - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. - real(pr), intent(out) :: half - real(pr) :: z_lim, cross_cn, z_cross, z_aux ! cn: carbon number - integer :: x_aux - - !A and B limit calculation. - z_lim = 0.0_pr - half = 0.506_pr ! modified by oscar, the variable half was removed fron the argument - !print*, scn - do while (z_lim.lt.plus_z) - half = half + 0.002d0 - cross_cn = scn(scn_nc)+half ! typically 19.508 in first try - z_cross = exp((a_blr*cross_cn) + b_blr) - a_lim = -(z_cross)/plus_z - b_lim = log(z_cross)-(a_lim*cross_cn) - x_aux=scn(scn_nc) - z_lim = 0._pr - do while (z_lim.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a_lim*x_aux+b_lim) - z_lim=z_lim+z_aux - end do - continue - end do - - c_max_lim = x_aux - !print*, a_lim,b_lim,c_max_lim - - end subroutine LimitLine - - subroutine Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) - !! This subroutine obtains the Cmax60 line constants for an oil - - implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the oil - integer, allocatable, intent(in) :: scn(:) !! set of singles cuts being considered in the oil - integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. - real(pr), intent(in) :: half - integer :: x_aux - real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range - real(pr) :: a_old, Zp_60, F, dF_dA - - cross_cn = scn(scn_nc)+half ! typically 19.51 - z_cross= exp((a_blr*cross_cn)+b_blr) - F_tol=1_pr - a_tol=1_pr - a_60=a_blr - var_range = 60.5_pr - cross_cn - - do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) - a_old=a_60 - Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 - F = plus_z - Zp_60 - dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected - a_60=a_old-(F/dF_dA) - a_tol=abs(a_old-a_60) - F_tol=abs(F) - end do - - b_60=log(z_cross)-a_60*cross_cn - z_sum=0.0_pr - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z) - x_aux=x_aux+1 - z_aux= exp(a_60*x_aux+b_60) - z_sum=z_sum+z_aux - end do - - c_max_60=x_aux - - !print*, a_60,b_60,c_max_60 - end subroutine Line_C60_max - - subroutine select_method(oil,mw_source,method,C,log_scn_z,plus_z,plus_mw) - - implicit none - type(FluidData) :: oil - !character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - real(pr), allocatable :: scn_mw(:) !! set to molecular weights calculated by \[M = 84-C(i-6)\] - real(pr) :: C !! C constants which is used in equation \[M = 84-C(i-6)\] - real(pr), allocatable :: scn_z(:) !! set of corresponding calculated mole fractions of scn cuts - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - real(pr), intent(out) :: plus_z !! composition of residual fraction from input file - real(pr) :: plus_mw !! calculated molecular weight of residual fraction - real(pr), allocatable, intent(out) :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cuts - real(pr), allocatable :: def_comp_moles(:) - real(pr), allocatable :: scn_moles(:) - real(pr) :: plus_moles - real(pr) :: total_moles - - - !oil = data_from_file(file) - - allocate(scn_z(oil%scn_nc)) - allocate(scn_mw(oil%scn_nc)) - allocate(log_scn_z(oil%scn_nc)) - allocate(def_comp_moles(oil%def_comp_nc)) - allocate(scn_moles(oil%scn_nc)) - - def_comp_moles = (oil%def_comp_w)/(oil%def_comp_mw) - - select case (mw_source) - - case("experimental") - scn_moles = (oil%scn_w)/(oil%scn_mw) - plus_moles = (oil%plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - log_scn_z = log(scn_z) - - case("calculated") - - if (method=="global_mw")then - scn_mw = 84 + C*(oil%scn-6) - scn_z = oil%product_z_mw_scn / scn_mw - sum_def_comp_z_plus = sum(oil%scn_z) + (oil%plus_z) - plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C - plus_mw = oil%product_z_mw_plus / plus_z - log_scn_z = log(scn_z) - end if - - if (method=="plus_mw")then - scn_mw = 84 + C*(oil%scn-6) - scn_moles = (oil%scn_w)/(scn_mw) - plus_moles = (oil%plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - log_scn_z = log(scn_z) - endif - - end select - - end subroutine select_method - - subroutine difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - - !! this subroutine ... - implicit none - type(FluidData) :: oil - logical :: start - !character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - integer :: c_max !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: scn_z(:) !! set of corresponding calculated mole fractions of scn cuts - integer :: n_init ! minimum carbon number obtained from the best linear regression - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cuts - real(pr) :: C !! C constants which is used in equation \[M = 84-C(i-6)\] - real(pr) :: plus_z !! composition of residual fraction from input file - real(pr) :: plus_mw !! calculated molecular weight of residual fraction - real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction - real(pr), intent(out) :: difference !! - real(pr) :: a !! output real variable. Slope of the best regression line. - real(pr) :: b !! output real variable. Intercept of the best regression line. - real(pr) :: r2 !! output real variable. Square correlation coefficient. - real(pr), allocatable :: x_blr(:), y_blr(:) - real(pr) :: a_blr, b_blr, a_60, b_60, a_lim, b_lim, half - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - real(pr), dimension(300) :: plus_z_i - real(pr), dimension(300) :: product_z_mw_plus_i - !real(pr) :: plus_mw - real(pr):: sum_z - integer :: j, k , k_old, n_best, x_aux, i_0! internal variables - integer :: c_max_blr, c_max_lim, c_max_60 - integer :: i - real(pr), allocatable :: scn_mw(:) !! set to molecular weights calculated by \[M = 84-C(i-6)\] - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - integer, dimension(300) :: carbon_number_plus - real(pr) :: denom - - !oil = data_from_file(file) - - allocate(scn_z(oil%scn_nc)) - allocate(scn_mw(oil%scn_nc)) - allocate(log_scn_z(oil%scn_nc)) - allocate(x_blr(oil%scn_nc)) - allocate(y_blr(oil%scn_nc)) - - 1 i_0 = oil%scn(oil%scn_nc) - - call select_method(oil,mw_source,method,C,log_scn_z,plus_z,plus_mw) ! select case - call Best_Linear_Regression(oil%scn_nc,oil%scn,log_scn_z,plus_z,a_blr,b_blr,r2, & - n_init,c_max_blr) - call LimitLine(oil%scn_nc,oil%scn,plus_z,a_blr,b_blr,a_lim,b_lim,c_max_lim,half) - call Line_C60_max (oil%scn_nc,oil%scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. - - if(a_blr < a_lim)then ! added by oscar 05/12/2023.se elimina por que se agregan las restriccones para C>12 y C<12. - a = a_lim - b = b_lim - else ! this line was removed call Line_C60_max - if(a_blr > a_60)then - a = a_60 - b = b_60 - else - a = a_blr - b = b_blr - end if - end if - - sum_z = 0.0_pr - i = 0 - do while (sum_z12 .and. C<14 )then !! mayor igual a 12 o menor igual 14. - denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) !! denominator of equation to obtain C value directly - C = (plus_z*(plus_mw-84))/denom - endif - - plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z - difference = plus_mw_cal-plus_mw - - end subroutine difference_mw_plus - - subroutine get_C_or_m_plus(oil,mw_source,method,start,C) - !! This subroutine... - implicit none - type(FluidData) :: oil - logical :: start - !character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - integer :: c_max !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cuts - real(pr) :: C !! C constants which is used in equation \[M = 84-C(i-6)\] - real(pr) :: plus_z !! composition of residual fraction from input file - real(pr) :: plus_mw !! calculated molecular weight of residual fraction - real(pr) :: difference !! - real(pr) :: difference_old, plus_mw_old - real(pr) :: C_old - real(pr) :: aux - - !oil = data_from_file(file) - - select case (mw_source) - - case("experimental") - C = 13.5 - Start = .true. - plus_mw = oil%plus_mw - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - print*, C, plus_mw - - case("calculated") - - if(method == 'global_mw') C=13 - if(method == 'plus_mw') C=14 - - Start = .true. - plus_mw = oil%plus_mw - - call difference_mw_plus(oil,mw_source,method,start,C,difference_old,plus_mw) - C_old = C - if(method == 'global_mw') C=min(12.8, C-0.07) - if(method == 'plus_mw') C=13.5 - Start = .false. - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - - do while (abs(difference) > 0.1) - aux = C - C = C - difference*(C-C_old)/(difference-difference_old) - C_old = aux - difference_old = difference - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - end do - print*, C, plus_mw - end select - - if(C>14.or.C<12)then - - if(C>14) C=14 - if(C<12) C=12 - - Start = .true. - plus_mw = oil%plus_mw - call difference_mw_plus(oil,mw_source,method,start,C,difference_old,plus_mw) - plus_mw_old = plus_mw - plus_mw = 0.9*plus_mw - Start = .false. - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - - do while (abs(difference)>0.00001) - aux = plus_mw - plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) - plus_mw_old = aux - difference_old = difference - call difference_mw_plus(oil,mw_source,method,start,C,difference,plus_mw) - end do - print*, C, plus_mw - end if - end subroutine get_C_or_m_plus - - end module pruebas - - - \ No newline at end of file diff --git a/src/routines.f90 b/src/routines.f90 deleted file mode 100644 index 445acec..0000000 --- a/src/routines.f90 +++ /dev/null @@ -1,768 +0,0 @@ -module routines - use constants - use dtypes, only: FluidData, FluidDataOut - use data_from_input, only: data_from_file, FluidData - -contains - - subroutine Linear_Regression(x, y, a, b, r2) - !! This subroutine computes the regression line for a data set of x, y variables. - implicit none - - integer :: i - integer :: n !! total number of data set. - real(pr), intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable - real(pr), intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable - real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables - real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line - real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line - real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient - - n = size(x) - !Calculation of a y b --> y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 - end do - - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 - end do - - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) - - end subroutine Linear_Regression - - subroutine Best_Linear_Regression(scn_nc, scn, scn_z, plus_z, a, b, r2, & - n_init, c_max_blr) - !! This subroutine calculates the best regression line for an fluid. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid - integer, intent(out) :: c_max_blr !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression - integer :: i, j, k , k_old, n_best, x_aux ! internal variables - real(pr), dimension(scn_nc) :: x_blr, y_blr ! x_blr: carbon number vector, y_blr: logarithm of mole fraction vector; for each step - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux ! auxiliary variables - - k=5 - r2=0.0001_pr - r2_old=0.00001_pr - r2_best=0.0001_pr - - do while (r2.gt.r2_old.or.r2_old.lt.0.9) - k_old=k - r2_old=r2 - a_old=a - b_old=b - - if (r2.gt.r2_best) then - r2_best=r2 - a_best=a - b_best=b - n_best=scn(scn_nc-k+2) - end if - - if (k.gt.scn_nc) then - if (r2.gt.r2_best)then - n_init=scn(scn_nc-k+2) - go to 22 - else - r2=r2_best - a=a_best - b=b_best - n_init=n_best - go to 22 - end if - end if - - j=1 - x_blr=0. - y_blr=0. - - do i=scn_nc-k+1, scn_nc - x_blr(j)=scn(i) - y_blr(j)=scn_z(i) - j=j+1 - end do - - call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) - k=k+1 - end do - - r2=r2_old - a=a_old - b=b_old - n_init=scn(scn_nc-k_old+2) - - 22 continue - - z_sum = 0d0 - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a*x_aux+b) - z_sum=z_sum+z_aux - end do - - c_max_blr=x_aux - - !print*, a, b, c_max_blr, r2, n_init - - end subroutine Best_Linear_Regression - - subroutine LimitLine(scn_nc, scn, plus_z, a_blr, b_blr, a_lim, b_lim, & - c_max_lim, half) - !! This subroutine obtains the limit line constants for an fluid. - implicit none - - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid. - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid. - integer, intent(out) :: c_max_lim !! maximum carbon number obtained for limit line. - !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr). - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression. - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression. - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file. - real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line. - real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. - real(pr), intent(out) :: half ! variable defined in the numerical method to converge the area under the curve of the function. - real(pr) :: z_lim, cross_cn, z_cross, z_aux - integer :: x_aux ! internal variable - - !A and B limit calculation. - z_lim = 0.0_pr - half = 0.506_pr ! modified by oscar, the variable half was removed from the argument - - do while (z_lim.lt.plus_z) - half = half + 0.002d0 - cross_cn = scn(scn_nc)+half ! typically 19.508 in first try - z_cross = exp((a_blr*cross_cn) + b_blr) - a_lim = -(z_cross)/plus_z - b_lim = log(z_cross)-(a_lim*cross_cn) - x_aux=scn(scn_nc) - z_lim = 0._pr - do while (z_lim.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a_lim*x_aux+b_lim) - z_lim=z_lim+z_aux - end do - continue - end do - - c_max_lim = x_aux - - end subroutine LimitLine - - subroutine Line_C60_max (scn_nc, scn, plus_z, a_blr, b_blr, half, a_60, & - b_60, c_max_60) - !! This subroutine obtains the Cmax60 line constants for an fluid. - - implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid - integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. - real(pr), intent(in) :: half !variable defined in the numerical method to converge the area under the curve of the function - integer :: x_aux ! auxiliary variable - real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range ! internal variables - real(pr) :: a_old, Zp_60, F, dF_dA ! internal variables - - cross_cn = scn(scn_nc)+half ! typically 19.51 - z_cross= exp((a_blr*cross_cn)+b_blr) - F_tol=1_pr - a_tol=1_pr - a_60=a_blr - var_range = 60.5_pr - cross_cn - - do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) - a_old=a_60 - Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 - F = plus_z - Zp_60 - dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected - a_60=a_old-(F/dF_dA) - a_tol=abs(a_old-a_60) - F_tol=abs(F) - end do - - b_60=log(z_cross)-a_60*cross_cn - z_sum=0.0_pr - x_aux=scn(scn_nc) - - do while (z_sum.lt.plus_z) - x_aux=x_aux+1 - z_aux= exp(a_60*x_aux+b_60) - z_sum=z_sum+z_aux - end do - - c_max_60=x_aux - - !print*, a_60,b_60,c_max_60 - end subroutine Line_C60_max - - subroutine select_method(fluid, mw_source, method, characterization) - !! This subroutine defines the calculation method to perform the characterization,& - !! based on the available experimental data. If the molecular weight data is experimental, - !! define the mole fractions, molecular weights and densities directly from the input file. - !! On the other hand, if the molecular weights are assumed, they are recalculated according - !! to the methodology described by Martin et al and the molar fractions of the fluid are recalculated. - - implicit none - type(FluidData), intent(inout) :: fluid !! Derived type to save inlet or experimental information of fluid to be characterized. - type(FluidDataOut), intent(inout) :: characterization !! Derive type variable to save results of characterization methodology. - character(len=*), intent(in), optional :: method - !! this option allows choose beetwen obtaining plus_mw or global_mw reported in the experimental information . - character(len=*), intent(in) :: mw_source !! This variable indicates the source of the input data, that is, whether they are experimental or not. - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - real(pr), allocatable :: def_comp_moles(:) !! - real(pr), allocatable :: scn_moles(:) !! - real(pr) :: plus_moles - real(pr) :: total_moles - - associate (& - scn_nc => fluid%scn_nc, def_comp_nc => fluid%def_comp_nc, & - def_comp_w => fluid%def_comp_w, def_comp_mw => fluid%def_comp_mw, & - scn_w => fluid%scn_w, plus_w => fluid%plus_w, scn => fluid%scn, & - plus_mw => characterization%plus_mw, scn_z => characterization%scn_z, & - plus_z => characterization%plus_z, log_scn_z => characterization%log_scn_z, & - C => characterization%C, scn_mw => characterization%scn_mw, & - scn_zm => characterization%scn_zm, & - plus_zm => characterization%plus_zm & - ) - - allocate(def_comp_moles(def_comp_nc)) - allocate(scn_moles(scn_nc)) - - def_comp_moles = (def_comp_w)/(def_comp_mw) - - select case (mw_source) - - case("experimental") - scn_mw = fluid%scn_mw - scn_moles = (scn_w)/(scn_mw) - plus_moles = (plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - - case("calculated") - - if (method=="global_mw")then - scn_mw = 84 + C*(scn-6) - scn_z = fluid%product_z_mw_scn / scn_mw - sum_def_comp_z_plus = sum(fluid%scn_z) + (fluid%plus_z) - plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C - plus_mw = fluid%product_z_mw_plus / plus_z - end if - - if (method=="plus_mw")then - scn_mw = 84 + C*(scn-6) - scn_moles = (scn_w)/(scn_mw) - plus_moles = (plus_w)/(plus_mw) - total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) - scn_z = scn_moles / total_moles - plus_z = plus_moles / total_moles - - endif - end select - scn_zm = scn_z*scn_mw - plus_zm = plus_z*plus_mw - log_scn_z = log(scn_z) - end associate - - end subroutine select_method - - subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & - characterization) - - !! this subroutine ... - implicit none - - type(FluidData), intent(inout) :: fluid - type(FluidDataOut), intent(inout) :: characterization - logical :: start - character(len=*), intent(in), optional :: method - character(len=*), intent(in) :: mw_source - real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction - real(pr), intent(out) :: difference - real(pr) :: half - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - real(pr):: sum_z - real(pr) :: denom - !! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] - real(pr) :: sum_def_comp_z_plus - !! compositions sum from (define component +1) to (plus component) - integer :: j, k , k_old, n_best, x_aux, i_0! internal variables - integer :: i - integer, dimension(300) :: carbon_number_plus - real(pr), dimension(300) :: plus_z_i - real(pr), dimension(300) :: product_z_mw_plus_i - real(pr), dimension(300) :: scn_i - - associate (& - scn_nc => fluid%scn_nc, scn => fluid%scn,& - plus_z => characterization%plus_z, plus_mw => characterization%plus_mw, & - log_scn_z => characterization%log_scn_z, C => characterization%C, & - a_blr => characterization%a_blr, b_blr => characterization%b_blr, & - a_60 => characterization%a_60, b_60 => characterization%b_60, & - a_lim => characterization%a_lim, b_lim => characterization%b_lim, & - r2 => characterization%r2, n_init => characterization%n_init, & - c_max_blr => characterization%c_max_lim, a => characterization%a, & - c_max_lim => characterization%c_max_lim, b => characterization%b, & - c_max_60 => characterization%c_max_lim, & - c_max => characterization%c_max & - ) - - 1 i_0 = scn(scn_nc) - - call select_method(fluid, mw_source, method, characterization) ! select case - call Best_Linear_Regression(scn_nc, scn, log_scn_z, & - plus_z, a_blr, b_blr, r2, n_init, c_max_blr) - call LimitLine(scn_nc, scn, plus_z, a_blr, b_blr,a_lim,b_lim,c_max_lim,half) - call Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. - - if(a_blr < a_lim)then - !Added by oscar 05/12/2023.se elimina por que se agregan las restricciones - !para characterization%C>12 y characterization%C<12. - a = a_lim - b = b_lim - else ! this line was removed call Line_C60_max - if(a_blr > a_60)then - a = a_60 - b = b_60 - else - a = a_blr - b = b_blr - end if - end if - - sum_z = 0.0_pr - i = 0 - do while (sum_z < plus_z .and. i< 300) - i = i+1 - plus_z_i(i) = exp(a*(i+i_0)+b) - sum_z = sum_z + plus_z_i(i) - product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) - carbon_number_plus(i) = i+i_0 - end do - - if (i == 300.and.sum_z < plus_z)then - if(start) C = C - 0.07 - if(.not.start) C = C - 0.01 - go to 1 - end if - - c_max = i+i_0 - plus_z_i(i) = plus_z_i(i) - (sum_z-plus_z) ! Adjustment to Zp (z20+) - product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) - - if (mw_source == "experimental" .and. C>12 .and. C<14 ) then - ! mayor igual a 12 o menor igual 14. - denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) - ! denominator of equation to obtain C value directly - C = (plus_z*(plus_mw-84))/denom - endif - - plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z - difference = plus_mw_cal-plus_mw - scn_i = [scn, carbon_number_plus(1:i)] - end associate - - characterization%plus_z_i = plus_z_i(1:i) - characterization%product_z_mw_plus_i = product_z_mw_plus_i(1:i) - characterization%carbon_number_plus = carbon_number_plus(1:i) - characterization%nc_plus = i - characterization%scn_i = scn_i - - - - end subroutine difference_mw_plus - - subroutine get_c_or_m_plus(fluid, mw_source, method, fix_C, characterization) - !! This subroutine... - implicit none - type(FluidData) :: fluid - type(FluidDataOut) :: characterization - logical :: start - character(len=*), intent(in) :: mw_source - character(len=*), intent(in), optional :: method - logical, intent(in) :: fix_C !! If C is < 12 or > 14 then fix to the - !! closest value. - !integer :: c_max !! output CN at which characterization%plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cutsn - real(pr) :: difference !! - real(pr) :: difference_old, plus_mw_old - real(pr) :: C_old - real(pr) :: aux - - associate (C => characterization%C, plus_mw => characterization%plus_mw) - - select case (mw_source) - case("experimental") - C = 13.5 - Start = .true. - plus_mw = fluid%plus_mw - call difference_mw_plus(fluid, mw_source, method, start, difference, & - characterization) - - case("calculated") - - if(method == 'global_mw') C=13 - if(method == 'plus_mw') C=14 - - Start = .true. - plus_mw = fluid%plus_mw - call difference_mw_plus(fluid, mw_source, method, start, & - difference_old, characterization) - C_old = C - if(method == 'global_mw') C = min(12.8, C-0.07) - if(method == 'plus_mw') C = 13.5 - Start = .false. - call difference_mw_plus(fluid, mw_source, method, start, difference, & - characterization) - - do while (abs(difference) > 0.1) - aux = C - C = C - difference*(C-C_old)/(difference-difference_old) - C_old = aux - difference_old = difference - call difference_mw_plus(fluid, mw_source, method, start, & - difference, characterization) - end do - end select - - if(fix_C .and. C>14 .or. C<12)then - - - if(C>14) C=14 - if(C<12) C=12 - - Start = .true. - plus_mw = fluid%plus_mw - call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & - difference_old, characterization) - plus_mw_old = plus_mw - plus_mw = 0.9*plus_mw - Start = .false. - call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & - difference, characterization) - - do while (abs(difference)>0.00001) - aux = plus_mw - plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) - plus_mw_old = aux - difference_old = difference - call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & - difference, characterization) - end do - end if - end associate - end subroutine get_c_or_m_plus - - subroutine density_funtion(fluid, mw_source, characterization) - !! this subroutine ... - implicit none - type(FluidData) :: fluid - type(FluidDataOut) :: characterization - character(len=*), intent(in) :: mw_source - !real(pr), allocatable :: volume_6plus_cal(:) - real(pr) :: volume_6plus_exp - real(pr) :: a_d_old, aux - real(pr) :: difference, difference_old - - - volume_6plus_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & - (fluid%product_z_mw_plus/fluid%plus_density) - - associate( & - a_d => characterization%a_d , & - b_d => characterization%b_d, & - volume_6plus_cal => characterization%volume_6plus_cal & - ) - ! Now find constants for the Density function - a_d = -0.50_pr ! initial guess - b_d = 0.685_pr - a_d*exp(-0.60_pr) - call calculate_volume_6plus(fluid, mw_source, characterization) - difference_old = volume_6plus_cal - volume_6plus_exp - a_d_old = a_d - a_d = -0.49_pr - b_d = 0.685_pr - a_d*exp(-0.60_pr) - call calculate_volume_6plus(fluid, mw_source, characterization) - difference = volume_6plus_cal - volume_6plus_exp - - do while (abs(difference) > 0.001_pr) - aux = a_d - a_d = a_d - difference*(a_d - a_d_old)/(difference - difference_old) - b_d = 0.685_pr - a_d*exp(-0.60_pr) - a_d_old = aux - difference_old = difference - call calculate_volume_6plus(fluid, mw_source, characterization) - difference = volume_6plus_cal - volume_6plus_exp - end do - - end associate - end subroutine density_funtion - - subroutine calculate_volume_6plus(fluid, mw_source, characterization) - !! this subroutine ... - implicit none - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - character(len=*), intent(in) :: mw_source - real(pr) :: volume_6plus_cal - real(pr), allocatable :: plus_density_cal(:), scn_density_cal(:), plus6_density(:) - - !allocate(density_plus_cal(0), density_scn_cal(0), plus6_density(0)) - - - associate(& - a_d => characterization%a_d, & - b_d => characterization%b_d, & - cn_plus => characterization%carbon_number_plus, & - z_m_plus_i => characterization%product_z_mw_plus_i, & - scn_z_i => characterization%scn_z, & - scn_mw_i => characterization%scn_mw & - ) - - scn_density_cal = ((a_d)*(exp(- (real(fluid%scn, pr))/10._pr))) + b_d - plus_density_cal =((a_d)*(exp(- (real(cn_plus, pr) )/10._pr))) + b_d - - - select case (mw_source) - case("experimental") - volume_6plus_cal = sum((scn_z_i*scn_mw_i)/(fluid%scn_density)) + & - sum(z_m_plus_i/plus_density_cal) - plus6_density = [fluid%scn_density, plus_density_cal] - - case("calculated") - volume_6plus_cal = sum((scn_z_i*scn_mw_i)/(scn_density_cal)) + & - sum(z_m_plus_i/plus_density_cal) - plus6_density = [scn_density_cal, plus_density_cal] - end select - - characterization%volume_6plus_cal = volume_6plus_cal - characterization%plus6_density = plus6_density - - end associate - - end subroutine calculate_volume_6plus - - subroutine lump (fluid, characterization) - !! This sobroutine ... - implicit none - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - integer :: last_C ! CN of last single cut: typically 19 - integer :: i_last ! Number of elements in the distribution of CN+ - integer :: last ! Total elements starting after defined components. - integer :: scn_nc_input !! inicial number of single cuts being considered in the oil from data input - integer :: scn_nc_new !! new number of single cuts being considered in the oil defined as from scn_nc_ps variable. - real(pr) :: plus_w_new - real(pr), allocatable :: z_m_plus_i(:) - real(pr), allocatable :: plus_z_i(:) - real(pr), allocatable :: var_aux_1(:) - real(pr), allocatable :: var_aux_2(:) - integer :: i, prev_i - integer, dimension(15) :: j_ps - integer :: i_ps ! auxiliary variables - real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps ! auxiliary variables - real(pr), allocatable :: plus_z_ps(:) - real(pr), allocatable :: plus_mw_ps(:) - integer :: numbers_ps - real(pr), allocatable :: density_ps(:) - real(pr), allocatable :: w_ps(:) - real(pr) :: sum_zm_last_ps - - associate (& - plus_zm => characterization%plus_zm, & - scn_zm => characterization%scn_zm, & - plus_z => characterization%plus_z, & - scn_z => characterization%scn_z, & - plus_w => fluid%plus_w, & - w => fluid%w, & - plus_mw => characterization%plus_mw, & - carbon_number_plus => characterization%carbon_number_plus, & - C => characterization%C & - ) - - last_C = fluid%scn(fluid%scn_nc) !19 - i_last = characterization%nc_plus !124 - last = fluid%scn_nc + i_last !138 - - allocate (z_m_plus_i(0)) - allocate (plus_z_i(0)) - allocate (var_aux_1(0)) - allocate (var_aux_2(0)) - allocate (plus_z_ps(0)) - allocate (plus_mw_ps(0)) - allocate(density_ps(0)) - allocate (w_ps(0)) - - scn_nc_input = fluid%scn_nc - scn_nc_new = fluid%scn_nc_ps - 6 - ! scn_nc_ps : CN from which all SCN fractions will be lumped - ! into the specified number of pseudos - - characterization%plus_w = plus_w - z_m_plus_i = characterization%product_z_mw_plus_i - plus_z_i = characterization%plus_z_i - - if (scn_nc_new < scn_nc_input) then - ! Plus Fraction needs to be extended to include lower CN's - plus_zm = plus_zm + sum(scn_zm(scn_nc_new + 1 : scn_nc_input)) - plus_z = plus_z + sum(scn_z(scn_nc_new + 1 : scn_nc_input)) !extended zp - plus_w_new = plus_w + sum(w(fluid%def_comp_nc+scn_nc_new + 1 & - : fluid%def_comp_nc+scn_nc_input)) - !! use new variable for plus_w because type fluidata can't be modified. - plus_mw = plus_zm / plus_z - z_m_plus_i = scn_zm(scn_nc_new + 1 : scn_nc_input ) - z_m_plus_i = [z_m_plus_i, characterization%product_z_mw_plus_i] - plus_z_i = scn_z(scn_nc_new + 1 : scn_nc_input ) - plus_z_i = [plus_z_i, characterization%plus_z_i] - i_last = i_last + scn_nc_input - scn_nc_new - last_C = last_C - scn_nc_input + scn_nc_new - characterization%plus_w = plus_w_new - characterization%product_z_mw_plus_i = z_m_plus_i - characterization%plus_z_i = plus_z_i - end if - - numbers_ps = fluid%numbers_ps - j_ps = 0 - ! Lumping into Nps pseudos - rec_zm = plus_zm / numbers_ps - ! Recommended value for the product z*M (proportional to weight) - !for each pseudo according to pedersen - remain_plus_zm = plus_zm - i_ps = 1._pr - sum_z = 0.0_pr - sum_zm = 0.0_pr - sum_volume_ps = 0.0_pr - - i = 0.0_pr - - do while (i_ps < numbers_ps) - i = i + 1 - j_ps(i_ps) = j_ps(i_ps) + 1 - sum_z = sum_z + plus_z_i(i) - sum_zm = sum_zm + z_m_plus_i(i) - sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / & - characterization%plus6_density(scn_nc_new+i) - if (z_m_plus_i(i+1) > 2*(rec_zm - sum_zm)) then - ! when adding one more would go too far - plus_z_ps = [plus_z_ps, sum_z] - plus_mw_ps = [plus_mw_ps, sum_zm/sum_z] - w_ps = [w_ps, (characterization%plus_w*(sum_zm/plus_zm))] - remain_plus_zm = remain_plus_zm - sum_zm - if (remain_plus_zm < rec_zm) numbers_ps = i_ps + 1 - carbon_number_plus(scn_nc_new + i_ps) = 6+((plus_mw_ps(i_ps)-84)/C) - density_ps = [density_ps, (sum_zm/sum_volume_ps)] - i_ps = i_ps + 1 - sum_z = 0._pr - sum_zm = 0._pr - sum_volume_ps = 0._pr - end if - end do - - plus_z_ps = [plus_z_ps, sum(plus_z_i(i+1:i_last))] - ! at this point, Nps is the order for the last NON-Asphaltenic pseudo comp. - ! (e.g. 4 if 5 is for Asphaltenes) - plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i+1:i_last))/plus_z_ps(numbers_ps)] - ! (zMp - sum(zMpi(1:i))) / zps(Nps) - w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps))* & - (plus_mw_ps(numbers_ps))/(plus_zm))] - carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & - ((plus_mw_ps(numbers_ps) - 84) / C) - sum_zm_last_ps = (sum(z_m_plus_i(i+1:i_last))) - density_ps = [density_ps, ((sum_zm_last_ps)/(sum(z_m_plus_i(i+1:i_last)) & - /(characterization%plus6_density(scn_nc_new+i+1:i_last))))] - j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) - - do i = 1, numbers_ps - print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) - end do - - end associate - - !characterization%product_z_mw_plus_i = z_m_plus_i - !characterization%plus_z_i = plus_z_i - characterization%last_C = last_C - characterization%i_last = i_last - - !! esto es lo que haria en el write, luego eliminar este comentario. - !prev_i = fluid%scn(1)-1 - !do i = 1, scn_nc_new - ! print*, prev_i + i , characterization%scn_z(i) , characterization%scn_mw(i) - !end do - - - end subroutine lump - - subroutine get_critical_constants(fluid, characterization) - !! this subroutine .... - type(FluidData), intent(in) :: fluid - type(FluidDataOut), intent(inout) :: characterization - !implicit none - - - - - - - - - end subroutine get_critical_constants - - type(FluidDataOut) function characterize(file, mw_source, method, fix_C) & - result(characterization) - - type(FluidData) :: fluid - character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method !! plus_mw or global_mw - character(len=*), intent(in) :: mw_source - logical, intent(in), optional :: fix_C - - - fluid = data_from_file(file) - allocate(characterization%scn_z(fluid%scn_nc)) - allocate(characterization%log_scn_z(fluid%scn_nc)) - allocate(characterization%scn_mw(fluid%scn_nc)) - allocate(characterization%carbon_number_plus(0)) - allocate(characterization%plus_z_i(0)) - allocate(characterization%product_z_mw_plus_i(0)) - allocate(characterization%scn_zm(fluid%scn_nc)) - allocate(characterization%scn_i(0)) - allocate(characterization%plus6_density(0)) - - - call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) - call density_funtion(fluid=fluid, mw_source=mw_source, characterization=characterization) - call lump(fluid=fluid, characterization=characterization) - - end function characterize - -end module routines - - diff --git a/src/tools/io.f90 b/src/tools/io.f90 new file mode 100644 index 0000000..7fde863 --- /dev/null +++ b/src/tools/io.f90 @@ -0,0 +1,22 @@ +module ftools__io + implicit none + + interface str + module procedure :: str_gen + end interface str + + contains + function str_gen(int_in) result(str_out) + class(*), intent(in) :: int_in + character(len=99) :: str_mid + character(len=:), allocatable :: str_out + + select type(int_in) + type is (real) + write (str_mid, *) int_in + type is (integer) + write (str_mid, *) int_in + end select + str_out = trim(adjustl(str_mid)) + end function str_gen + end module ftools__io \ No newline at end of file diff --git a/src/types.f90 b/src/types.f90 index afec0a3..242f369 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -7,30 +7,46 @@ module dtypes integer :: scn_nc !! number of single cuts being considered in the oil integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos integer :: numbers_ps !! number of pseudos in which the scn fractions grouped. - character(len=:), allocatable :: filename - integer, allocatable :: scn (:) !! set of singles cuts being considered in the oil - character(len=15), allocatable :: def_components (:) !! set of defined components being considered in the oil - character(len=15) :: scn_plus !! name of residual fraction - real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components - real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), allocatable :: def_comp_mw(:) !! set of corresponding molecular weights of defined components + logical :: density_setup + integer :: number_plus_density + character(len=:), allocatable :: filename + integer, allocatable :: scn (:) !! set of singles cuts being considered in the oil + character(len=15), allocatable :: def_components (:) !! set of defined components being considered in the oil + character(len=15) :: scn_plus !! name of residual fraction + real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components + real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts + real(pr), allocatable :: def_comp_mw(:) !! set of corresponding molecular weights of defined components real(pr), allocatable :: scn_mw(:) !! set of corresponding molecular weights of scn cuts real(pr), allocatable:: product_z_mw_def_comp(:) !! product between composition and molecular weight of defind components real(pr), allocatable:: product_z_mw_scn(:) !! product between composition and molecular weight of scn fractions real(pr) :: sum_z_mw_i !! sum of the product between composition and molecular weight of the fluid's compounds real(pr), allocatable :: w(:) !! mass fractions of the fluid's compounds real(pr) :: plus_z !! composition of residual fraction - real(pr) :: plus_mw !! molecular weight of residual fraction + real(pr) :: plus6_z_exp !! composition of residual fraction + real(pr) :: plus7_z_exp !! composition of residual fraction + real(pr) :: plus12_z_exp !! composition of residual fraction + real(pr) :: plus30_z_exp !! composition of residual fraction + real(pr) :: plus_mw !! molecular weight of residual fraction + real(pr) :: plus6_mw_exp !! molecular weight of residual fraction + real(pr) :: plus7_mw_exp !! molecular weight of residual fraction + real(pr) :: plus12_mw_exp !! molecular weight of residual fraction + real(pr) :: plus30_mw_exp !! molecular weight of residual fraction real(pr) :: product_z_mw_plus !! product between composition and molecular weight of residual fraction real(pr), allocatable :: def_comp_w(:) !! mass fractions of the defined compounds real(pr), allocatable :: scn_w(:) !! !! mass fractions of the scn-s compounds real(pr) :: plus_w !! mass fractions of the plus fraction real(pr), allocatable :: scn_density(:) !! set of corresponding densities of scn cuts real(pr) :: plus_density !! experimental density of the plus fraction + real(pr) :: plus6_density_exp + real(pr) :: plus7_density_exp + real(pr) :: plus12_density_exp + real(pr) :: plus30_density_exp + end type FluidData type :: FluidDataOut ! incluir aqui las variables que quiero que salgan como salida + type(FluidData) :: input_data real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts calculated real(pr), allocatable :: log_scn_z(:) !! set of corresponding mole fractions of scn cuts calculated real(pr) :: plus_mw !! molecular weight of residual fraction @@ -57,14 +73,114 @@ module dtypes real(pr), allocatable :: product_z_mw_plus_i(:) real(pr) :: a_d !! ad constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. real(pr) :: b_d !! bd constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. - real(pr) :: volume_6plus_cal + real(pr) :: volume_cal + real(pr) :: volume_exp integer :: last_C integer :: i_last + integer :: last + integer :: scn_nc_new real(pr), allocatable :: scn_i(:) real(pr), allocatable :: plus6_density(:) real(pr) :: plus_w + real(pr) :: plus_density + real(pr), allocatable :: mol_fraction(:) + real(pr), allocatable :: lumped_z(:) + real(pr), allocatable :: lumped_mw(:) + real(pr), allocatable :: lumped_densities(:) + real(pr), allocatable :: critical_temperature(:) + real(pr), allocatable :: critical_pressure(:) + real(pr), allocatable :: acentric_factor(:) + real(pr), allocatable :: m_funtion (:) + + contains + private + procedure, pass :: write => write_result + generic, public :: write (FORMATTED) => write end type FluidDataOut +contains + + subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) + use ftools__io, only: str + use defined_critical_parameters + + implicit none + class(FluidDataOut), intent(in) :: characterization + integer :: i, i_prev + integer, intent(in) :: unit + integer, intent(out) :: iostat + character(*), optional, intent(in) :: iotype + character(*), optional, intent(inout) :: iomsg + integer, optional, intent(in) :: v_list(:) + character(len=*), parameter :: str_fmt_1 = "(A4,7(A15,2x),/)" + character(len=*), parameter :: num_fmt_1 = "(A3,1x,5(E15.5,2x),/)" + character(len=*), parameter :: num_fmt_2 = "(A3,1x,7(E15.5,2x),/)" + character(len=*), parameter :: str_fmt_2 = "(130('-'),/),/)" + + associate(& + def_nc => characterization%input_data%def_comp_nc, & + def_name => characterization%input_data%def_components, & + def_mw => characterization%input_data%def_comp_mw, & + z => characterization%mol_fraction, & + mw => characterization%lumped_mw, & + rho => characterization%lumped_densities, & + scn => characterization%input_data%scn, & + scn_nc_new => characterization%scn_nc_new, & + tc => characterization%critical_temperature, & + pc => characterization%critical_pressure, & + omega => characterization%acentric_factor, & + m_funtion => characterization%m_funtion & + ) + + + + write(unit, fmt=str_fmt_2, iostat = iostat) + !write(*,"(A64,/)") "----------------------------------------------------------------" + !write(*,*) "Best feasible regresion parameters" + !write(unit, fmt=str_fmt_1, iostat = iostat) "" + !write(unit, fmt=num_fmt_2, iostat = iostat) "Init_BFR:", characterization%n_init + !write(*,*) "A:", characterization%a," ","B:", characterization%b + !write(*,*) "C:", characterization%C + !write(*,*) "MW+:", characterization%plus_mw + !write(*,*) "Cmax:", characterization%c_max + !write(*,*) "" + !write(*,*) "Density funtion parameters parameters" + !write(*,*) "" + !write(*,*) "ad:", characterization%a_d," ","bd:", characterization%b_d + !write(*,*) "----------------------------------------------------------------" + !write(*,*) "" + !write(*,*) "compositional result" + !write(*,*) "" + + + + + + + + write(unit, fmt=str_fmt_1, iostat = iostat) "Comp", "Z", "Mw", "Tc", & + "Pc", "Omega", "Rho", "M" + ! defined components + do i = 1, def_nc + write(unit, fmt=num_fmt_1, iostat = iostat) def_name(i), z(i), & + def_mw(i), tc_def(i), pc_def(i), om_def(i) + end do + ! + i_prev = scn(1) - 1 + do i = 1, scn_nc_new + write(unit, fmt=num_fmt_2, iostat = iostat) adjustl('C' // str(i_prev + i)), & + z(i+def_nc), mw(i), tc(i), pc(i), omega(i), rho(i), m_funtion(i) + end do + ! + do i = 1 + scn_nc_new, size(mw) + write(unit, fmt=num_fmt_2, iostat = iostat) 'ps' // str(i - scn_nc_new), & + z(i+def_nc), mw(i), rho(i), tc(i), pc(i), omega(i), m_funtion(i) + end do + + end associate + end subroutine write_result + + end module dtypes