Merge pull request #691 from bobmyhill/refactor_solver

refactor nonlinear solver

Author: bobmyhill

burnman/classes/solutionmodel.py

@@ -6,14 +6,13 @@
 
 import importlib
 import numpy as np
+from ..constants import logish_eps
 from ..utils.chemistry import process_solution_chemistry
 from .. import constants
 import warnings
 from .material import material_property, cached_property
 from dataclasses import make_dataclass
 
-np_eps = np.finfo(float).eps
-
 Interaction = make_dataclass(
     "Interaction", ["inds", "expts", "f_r", "m_jr", "f_rs", "m_jrs"]
 )
@@ -106,26 +105,26 @@ def _non_ideal_interactions_subreg(p, n_endmembers, Wijk):
     return Wint
 
 
-def logish(x, eps=np_eps):
+def logish(x, eps=logish_eps):
     """
     2nd order series expansion of log(x) about eps:
     log(eps) - sum_k=1^infty (f_eps)^k / k
     Prevents infinities at x=0
     """
-    f_eps = 1.0 - x / eps
-    mask = x > eps
-    ln = np.where(x <= eps, np.log(eps) - f_eps - f_eps * f_eps / 2.0, 0.0)
+    f_eps = 1.0 - x / eps[0]
+    mask = x > eps[0]
+    ln = np.where(x <= eps[0], np.log(eps[0]) - f_eps - f_eps * f_eps / 2.0, 0.0)
     ln[mask] = np.log(x[mask])
     return ln
 
 
-def inverseish(x, eps=np_eps):
+def inverseish(x, eps=logish_eps):
     """
     1st order series expansion of 1/x about eps: 2/eps - x/eps/eps
     Prevents infinities at x=0
     """
-    mask = x > eps
-    oneoverx = np.where(x <= eps, 2.0 / eps - x / eps / eps, 0.0)
+    mask = x > eps[0]
+    oneoverx = np.where(x <= eps[0], 2.0 / eps[0] - x / eps[0] / eps[0], 0.0)
     oneoverx[mask] = 1.0 / x[mask]
     return oneoverx
 

burnman/constants.py

@@ -3,6 +3,7 @@
 # GPL v2 or later.
 
 import scipy.constants
+import numpy as np
 
 """
 molar gas constant (R) in J mol^-1 K^-1
@@ -42,3 +43,10 @@
 1 cm^-1 in J/mol
 """
 invcm = 11.9627
+
+
+"""
+Small number for numerical stability
+Kept in list form for mutability
+"""
+logish_eps = [np.finfo(float).eps]