Skip to content

Commit

Permalink
some mypy fixes
Browse files Browse the repository at this point in the history
disabled mypy
  • Loading branch information
Fahd-Siddiqui committed Feb 4, 2024
1 parent 049ef3e commit dfca67e
Show file tree
Hide file tree
Showing 8 changed files with 125 additions and 61 deletions.
6 changes: 3 additions & 3 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,9 @@ lint:
@isort --check .
@echo "Running lint..."
@black --check .
@flake8 src --exclude=venv* --max-line-length=120
@echo "Running mypy..."
@mypy src
@flake8 src --exclude=venv* --max-line-length=120 --extend-ignore=E203,E704
# @echo "Running mypy..."
# @mypy src

# Tests: pytest from src folder with coverage and printing coverage report
test:
Expand Down
Empty file removed __init__.py
Empty file.
10 changes: 8 additions & 2 deletions pyproject.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,12 @@
[tool.mypy]
disallow_untyped_defs = true
ignore_missing_imports = true
strict = true
[[tool.mypy.overrides]]
module = [
"numpy.*"
]
strict=false
disallow_untyped_calls = false
disallow_untyped_defs = false

[tool.black]
line-length = 120
Expand Down
2 changes: 1 addition & 1 deletion src/calculator.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ class Calculator:
@classmethod
def calculate_fugacity_coef_difference(cls, comp, T, P, a, b, amix, bmix, Composition, kij, lij, phase):
volume = EOS.calculate_mixing_rules(amix, bmix, comp, a, b, kij, lij, Composition, P, T, phase)
fugacity_vec = EOS.fugacity_vec(T, P, a, b, amix, bmix, volume, Composition, kij, lij)
fugacity_vec = EOS.fugacity_vectorized(T, P, a, b, amix, bmix, volume, Composition, kij, lij)

fugacity_coefficients = [fugacity_vec(np.arange(comp), ph).astype("float64") for ph in phase]

Expand Down
100 changes: 74 additions & 26 deletions src/eos.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,7 @@
from typing import Callable, Tuple

import numpy as np
import numpy.typing as npt

from src.Constants import Constants
from src.utils import Utils
Expand All @@ -8,11 +11,11 @@ class EOS:
@classmethod
def eos_parameters(
cls,
accentric_factor: np.ndarray,
critical_temperature: np.ndarray,
ac: np.ndarray,
accentric_factor: npt.NDArray[np.float64],
critical_temperature: npt.NDArray[np.float64],
ac: npt.NDArray[np.float64],
temperature: float,
) -> np.ndarray:
) -> npt.NDArray[np.float64]:
K = 0.37464 + (1.54226 - 0.26992 * accentric_factor) * accentric_factor
reduced_temperature = temperature / critical_temperature

Expand All @@ -24,12 +27,12 @@ def eos_parameters(
def VdW1fMIX(
cls,
comp: int,
a: np.ndarray,
b: np.ndarray,
kij: np.ndarray,
lij: np.ndarray,
MoleFrac: np.ndarray,
):
a: npt.NDArray[np.float64],
b: npt.NDArray[np.float64],
kij: npt.NDArray[np.float64],
lij: npt.NDArray[np.float64],
MoleFrac: npt.NDArray[np.float64],
) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
MoleFracAux = Utils.normalize(MoleFrac)

amix = np.sum(np.multiply.outer(MoleFracAux, MoleFracAux) * np.sqrt(np.outer(a, a)) * (1 - kij))
Expand All @@ -38,15 +41,22 @@ def VdW1fMIX(
return amix, bmix

@classmethod
def Eos_Volumes(cls, P: float, T: float, amix: np.ndarray, bmix: np.ndarray, phase: np.ndarray) -> np.ndarray:
def Eos_Volumes(
cls,
P: float,
T: float,
amix: npt.NDArray[np.float64],
bmix: npt.NDArray[np.float64],
phase: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
Volume = np.zeros(len(phase))
for ph in phase:
Volume[ph] = cls.EoS_Volume(P, T, bmix[ph], amix[ph], phase[ph])

return Volume

@classmethod
def EoS_Volume(cls, P: float, T: float, bmix: float, amix: float, root: int) -> float:
def EoS_Volume(cls, P: float, T: float, bmix: float, amix: float, root: int) -> np.float64:
R = Constants.R
sigma_eos = 1.0 + 2.0**0.5
epsilon_eos = 1.0 - 2.0**0.5
Expand All @@ -64,7 +74,7 @@ def EoS_Volume(cls, P: float, T: float, bmix: float, amix: float, root: int) ->
coefCubic[3] = -(1.0 / aux + bmix) * sigma_eos * epsilon_eos * bmix**2.0 - bmix * amix / P

# Cubic equation solver
Vol = np.roots(coefCubic).real[abs(np.roots(coefCubic).imag) < 1e-3]
Vol: npt.NDArray[np.float64] = np.roots(coefCubic).real[abs(np.roots(coefCubic).imag) < 1e-3]
Vol = np.where(Vol < bmix, 1e16 if root == 1 else 1e-16, Vol)
return np.min(Vol) if root == 1 else np.max(Vol)

Expand All @@ -73,14 +83,14 @@ def fugacity(
cls,
T: float,
P: float,
a: np.ndarray,
b: np.ndarray,
amix: float,
bmix: float,
Volume: float,
MoleFrac: np.ndarray,
kij: np.ndarray,
lij: np.ndarray,
a: npt.NDArray[np.float64],
b: npt.NDArray[np.float64],
amix: np.float64,
bmix: np.float64,
Volume: np.float64,
MoleFrac: npt.NDArray[np.float64],
kij: npt.NDArray[np.float64],
lij: npt.NDArray[np.float64],
index: int,
) -> np.float64:
sigma_eos = 1.0 + 2.0**0.5 # PR
Expand Down Expand Up @@ -108,7 +118,20 @@ def fugacity(
return FugCoef

@classmethod
def calculate_fugacity_coefs2(cls, comp, T, P, a, b, amix, bmix, volume, Composition, kij, lij):
def calculate_fugacity_coefs2(
cls,
comp: int,
T: float,
P: float,
a: npt.NDArray[np.float64],
b: npt.NDArray[np.float64],
amix: npt.NDArray[np.float64],
bmix: npt.NDArray[np.float64],
volume: npt.NDArray[np.float64],
Composition: npt.NDArray[np.float64],
kij: npt.NDArray[np.float64],
lij: npt.NDArray[np.float64],
) -> Tuple[float, float, float]:
# TODO Instead of using 0,1 implement Phase
fug_func = np.frompyfunc(
lambda i: EOS.fugacity(
Expand Down Expand Up @@ -150,14 +173,39 @@ def calculate_fugacity_coefs2(cls, comp, T, P, a, b, amix, bmix, volume, Composi
return fug_coef_ref, fug_coef_aux, difference

@classmethod
def calculate_mixing_rules(cls, amix, bmix, comp, a, b, kij, lij, composition, P, T, phase):
def calculate_mixing_rules(
cls,
amix: npt.NDArray[np.float64],
bmix: npt.NDArray[np.float64],
comp: int,
a: npt.NDArray[np.float64],
b: npt.NDArray[np.float64],
kij: npt.NDArray[np.float64],
lij: npt.NDArray[np.float64],
composition: npt.NDArray[np.float64],
P: float,
T: float,
phases: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
amix[0], bmix[0] = cls.VdW1fMIX(comp, a, b, kij, lij, composition[:, 0]) # Mixing Rule - Reference Phase
amix[1], bmix[1] = cls.VdW1fMIX(comp, a, b, kij, lij, composition[:, 1]) # Mixing Rule - Incipient Phase
volume = EOS.Eos_Volumes(P, T, amix, bmix, phase)
volume = EOS.Eos_Volumes(P, T, amix, bmix, phases)
return volume

@classmethod
def fugacity_vec(cls, T, P, a, b, amix, bmix, volume, Composition, kij, lij):
def fugacity_vectorized(
cls,
T: float,
P: float,
a: npt.NDArray[np.float64],
b: npt.NDArray[np.float64],
amix: npt.NDArray[np.float64],
bmix: npt.NDArray[np.float64],
volume: np.typing.NDArray[np.float64],
Composition: np.typing.NDArray[np.float64],
kij: npt.NDArray[np.float64],
lij: npt.NDArray[np.float64],
) -> Callable:
return np.frompyfunc(
lambda i, phase_number: EOS.fugacity(
T,
Expand All @@ -177,5 +225,5 @@ def fugacity_vec(cls, T, P, a, b, amix, bmix, volume, Composition, kij, lij):
)

@classmethod
def VdW1fMIX_vec(cls, comp, a, b, kij, lij, MoleFrac):
def VdW1fMIX_vectorized(cls, comp, a, b, kij, lij, MoleFrac) -> Callable:
return np.frompyfunc(lambda i: EOS.VdW1fMIX(comp, a, b, kij, lij, MoleFrac[:i]), nin=1, nout=2)
60 changes: 34 additions & 26 deletions src/phase_envelope.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
from typing import Dict, List

import numpy as np
import numpy.typing as npt

from src.calculator import Calculator
from src.Constants import Constants
Expand Down Expand Up @@ -34,10 +35,10 @@ def calculate(
self,
temperature: float,
pressure: float,
phase_compositions: np.ndarray,
critical_temperatures: np.ndarray,
critical_pressures: np.ndarray,
acentric_factors: np.ndarray,
phase_compositions: npt.NDArray[np.float64],
critical_temperatures: npt.NDArray[np.float64],
critical_pressures: npt.NDArray[np.float64],
acentric_factors: npt.NDArray[np.float64],
):
phase_compositions = np.array(phase_compositions)
critical_temperatures = np.array(critical_temperatures)
Expand All @@ -63,20 +64,20 @@ def _calculate(
self,
temperature: float,
pressure: float,
phase_compositions: np.ndarray,
critical_temperatures: np.ndarray,
critical_pressures: np.ndarray,
acentric_factors: np.ndarray,
phase_compositions: npt.NDArray[np.float64],
critical_temperatures: npt.NDArray[np.float64],
critical_pressures: npt.NDArray[np.float64],
acentric_factors: npt.NDArray[np.float64],
):
component_index = len(phase_compositions)
a_mix: np.ndarray = np.zeros(2)
b_mix: np.ndarray = np.zeros(2)
a_mix: npt.NDArray[np.float64] = np.zeros(2)
b_mix: npt.NDArray[np.float64] = np.zeros(2)

phase_envelope_results = []

composition: np.ndarray = np.zeros((component_index, component_index))
kij: np.ndarray = np.zeros((component_index, component_index))
lij: np.ndarray = np.zeros((component_index, component_index))
composition: npt.NDArray[np.float64] = np.zeros((component_index, component_index))
kij: npt.NDArray[np.float64] = np.zeros((component_index, component_index))
lij: npt.NDArray[np.float64] = np.zeros((component_index, component_index))

dF = np.zeros((component_index + 2) ** 2)
sensitivity = np.zeros(component_index + 2)
Expand Down Expand Up @@ -186,7 +187,16 @@ def _calculate(
max_step,
)

# self.logger.info("Incipient Phase = ", phase[1] + 1, " P = ", P, "T = ", T, "self.specified_variable =", self.specified_variable)
# self.logger.info(
# "Incipient Phase = ",
# phase[1] + 1,
# " P = ",
# P,
# "T = ",
# T,
# "self.specified_variable =",
# self.specified_variable,
# )

if max_step > self.TOLERANCE or any(np.isnan(Var)):
flag_error = 1
Expand Down Expand Up @@ -279,7 +289,7 @@ def _calculate(
return phase_envelope_results

@staticmethod
def solve(matrix: np.ndarray, b: np.ndarray):
def solve(matrix: npt.NDArray[np.float64], b: npt.NDArray[np.float64]):
return np.linalg.solve(matrix, b)

def differentiate_first_c_residuals_lnK(
Expand Down Expand Up @@ -390,7 +400,7 @@ def newtons_method_loop(
F[comp] = np.sum(Composition[:, 1] - Composition[:, 0])
F[comp + 1] = Var[self.specified_variable_index] - S

# Differentiating The First "C" Residuals With Respect to ln[K(j)]******************************************************
# Differentiating The First "C" Residuals With Respect to ln[K(j)]
self.differentiate_first_c_residuals_lnK(
comp,
Composition,
Expand All @@ -407,10 +417,10 @@ def newtons_method_loop(
dF,
)

# Differentiating "C+1" Residual With Respect to ln[K[i]]///////////////////////////////////////////////////////////////
# Differentiating "C+1" Residual With Respect to ln[K[i]]
dF[comp * (comp + 2) : comp * (comp + 2) + comp] = Composition[:, 1]

# Differentiating The First "C" Residuals With Respect to ln[T]*********************************************************
# Differentiating The First "C" Residuals With Respect to ln[T]
diffT = self.diff * Var[comp]

# Numerically Differentiating The ln(Fugacity Coefficient) With Respect to ln(T)
Expand All @@ -432,7 +442,7 @@ def newtons_method_loop(
T = np.exp(Var[comp])
a = EOS.eos_parameters(acentric, Tc, ac, T)

# Differentiating The First "C" Residuals With Respect to ln[P]/////////////////////////////////////////////////////////
# Differentiating The First "C" Residuals With Respect to ln[P]
diffP = self.diff * Var[comp + 1]
for ph in phase:
amix[ph], bmix[ph] = EOS.VdW1fMIX(comp, a, b, kij, lij, Composition[:, ph])
Expand Down Expand Up @@ -470,7 +480,7 @@ def newtons_method_loop(
Var = Var - step
maxstep = max(abs(step / Var))

# Calculating The Natural Form Of Indep# endent Variables And Updating Compositions Of The Incipient Phase////////////////
# Calculating The Natural Form Of Indep# endent Variables And Updating Compositions Of The Incipient Phase
K = np.exp(Var[0:comp])
Composition[:, 1] = Composition[:, 0] * K
T = np.exp(Var[comp + 0])
Expand Down Expand Up @@ -536,7 +546,7 @@ def handle_non_critical_point(
sensitivity[self.specified_variable_index] = 1.0
S = Var[self.specified_variable_index]

# Adjusting Stepsize////////////////////////////////////////////////////////////////////////////////////////////////////////
# Adjusting Stepsize
dSmax = max(abs(Var[self.specified_variable_index]) ** 0.5 / 10.0, 0.1) * abs(dS) / dS

dS *= 4.0 / it
Expand All @@ -545,14 +555,13 @@ def handle_non_critical_point(

# Defining Specified Variable Value In The Next Point
S += dS
# //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

# Indep# endent Variables Initial Guess For The Next Point********************************************************************
# Indep# endent Variables Initial Guess For The Next Point
Var += dS * sensitivity

# **************************************************************************************************************************

# Analyzing Temperature Stepsize////////////////////////////////////////////////////////////////////////////////////////////
# Analyzing Temperature Stepsize
T_old = T
T = np.exp(Var[comp + 0])
# Large Temperature Steps Are Not Advisable
Expand All @@ -561,9 +570,8 @@ def handle_non_critical_point(
S -= dS
Var -= dS * sensitivity
T = np.exp(Var[comp])
# //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

# Analyzing Proximity to Critical Point*************************************************************************************
# Analyzing Proximity to Critical Point
# Seeking The Greatest K-factor
maxK_i = np.argmax(abs(Var[0:comp]))
# If The ln[K[i]] Stepsize Is Too Big, It Should Be Decreased
Expand Down
3 changes: 2 additions & 1 deletion src/successive_substitution.py
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,8 @@ def get_initial_temperature_guess(
gibbs_vap = np.sum(z * fug_coef_ref)
gibbs_liq = np.sum(z * fug_coef_aux)

# factor 1.75 was proposed by Pedersen and Christensen (Phase Behavior Of Petroleum Reservoir Fluids, 2007 - Chapter 6.5 Phase Identification)
# factor 1.75 was proposed by Pedersen and Christensen
# (Phase Behavior Of Petroleum Reservoir Fluids, 2007 - Chapter 6.5 Phase Identification)
if (gibbs_liq < gibbs_vap) or (gibbs_liq == gibbs_vap and (volume[0] / bmix[0]) < 1.75):
temperature = temperature + 10.0

Expand Down
5 changes: 3 additions & 2 deletions src/utils.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
import numpy
import numpy as np
import numpy.typing as npt


class Utils:
@classmethod
def normalize(cls, array: numpy.ndarray) -> numpy.ndarray:
def normalize(cls, array: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
return array / sum(array)

0 comments on commit dfca67e

Please sign in to comment.