diff --git a/burnman/classes/solutionmodel.py b/burnman/classes/solutionmodel.py index 3761c276..a19ce2fb 100644 --- a/burnman/classes/solutionmodel.py +++ b/burnman/classes/solutionmodel.py @@ -13,7 +13,12 @@ import sparse import string from copy import deepcopy -from .material import material_property +from .material import material_property, cached_property +from dataclasses import make_dataclass + +Interaction = make_dataclass( + "Interaction", ["inds", "expts", "f_r", "m_jr", "f_rs", "m_jrs"] +) try: ag = importlib.import_module("autograd") @@ -127,6 +132,41 @@ def inverseish(x, eps=1.0e-5): return oneoverx +def dpdn(molar_amounts, n, ones, eye): + """ + The partial derivative of endmember proportions + with respect to endmember amounts. + + :param molar_amounts: molar amounts of independent endmembers + :type molar_fractions: 1D numpy array + :param n: sum of the endmember amounts (usually equal to one) + :type n: float + :param ones: a vector of ones of length equal to the number of endmembers + :type ones: 1D numpy array + :param eye: the identity matrix of size equal to the number of endmembers + :type eye: 2D numpy array + """ + return (eye - np.einsum("k, m->km", molar_amounts, ones)) / n + + +def d2pdndn(molar_amounts, nsqr, ones, eyeones): + """ + The second partial derivative of endmember proportions + with respect to endmember amounts. + + :param molar_amounts: molar amounts of independent endmembers + :type molar_fractions: 1D numpy array + :param nsqr: square of the sum of the endmember amounts (usually equal to one) + :type n: float + :param ones: a vector of ones of length equal to the number of endmembers + :type ones: 1D numpy array + :param eyeones: delta_ij 1_k + delta_ik 1_j + :type eyeones: 3D numpy array + """ + + return (2.0 * np.einsum("n, m, k->knm", ones, ones, molar_amounts) - eyeones) / nsqr + + class SolutionModel(object): """ This is the base class for a solution model, intended for use @@ -506,6 +546,37 @@ def activity_coefficients(self, pressure, temperature, molar_fractions): def activities(self, pressure, temperature, molar_fractions): return self._ideal_activities(molar_fractions) + @cached_property + def ones(self): + """ + A vector of ones with length equal to the number of endmembers + :return: ones + :rtype: 1D numpy array + """ + return np.ones(self.n_endmembers) + + @cached_property + def eye(self): + """ + An identity matrix with size equal to the number of endmembers + :return: eye + :rtype: 2D numpy array + """ + return np.eye(self.n_endmembers) + + @cached_property + def eyeones(self): + """ + A convenience function consisting of two concatenations + of an identity matrix and ones vector with size + equal to the number of endmembers. + :return: delta_ij 1_k + delta_ik 1_j + :rtype: 3D numpy array + """ + return np.einsum("km, n->kmn", self.eye, self.ones) + np.einsum( + "kn, m->kmn", self.eye, self.ones + ) + class AsymmetricRegularSolution(IdealSolution): """ @@ -1025,9 +1096,9 @@ class PolynomialSolution(IdealSolution): :param ESV_interactions: A list containing lists where the first three elements are energy, entropy and volume interactions and the rest of the elements are indices of the transformed endmembers to which those - interactions correspond. - For example, [2., 0., 0., 0, 1, 1] would correspond to an interaction - of 2*p'[0]*p'[1]*p'[1]. + interactions correspond, and the proportion exponents. + For example, [2., 0., 0., 0, 4, 1, 1] would correspond to an interaction + of 2*p'[0]^4*p'[1]^1. :type ESV_interactions: list of lists :param interaction_endmembers: A list of minerals involved in the interaction terms. :type interaction_endmembers: list of :class:`burnman.Mineral` objects @@ -1036,8 +1107,8 @@ class PolynomialSolution(IdealSolution): coefficients for each of the interaction_endmembers and the rest of the elements are indices of the transformed endmembers to which those interactions correspond. - For example, [1., 0., -1., 0, 1, 1] would correspond to an interaction - of (mbr[0].gibbs - mbr[2].gibbs)*p'[0]*p'[1]*p'[1]. + For example, [1., 0., -1., 0, 4, 1, 1] would correspond to an interaction + of (mbr[0].gibbs - mbr[2].gibbs)*p'[0]^4*p'[1]^1. :type endmember_coefficients_and_interactions: list of lists :param transformation_matrix: The interactions for a given solution may be most compactly expressed not as a polynomial function of the @@ -1062,19 +1133,29 @@ def __init__( self.endmembers = endmembers self.W_ESV = None + self.c_ESV = None + self.n_W_ESV = 0 if ESV_interactions is not None: - self.W_ESV = self.make_interaction_arrays(ESV_interactions) + self.W_ESV, self.c_ESV = self._make_interaction_arrays(ESV_interactions, 3) + self.n_W_ESV = len(ESV_interactions) self.n_interaction_endmembers = len(interaction_endmembers) self.interaction_endmembers = interaction_endmembers self.W_mbr = None + self.c_mbr = None + self.n_W_mbr = 0 if self.n_interaction_endmembers > 0: - self.W_mbr = self.make_endmember_interaction_arrays( - endmember_coefficients_and_interactions + self.W_mbr, self.c_mbr = self._make_interaction_arrays( + endmember_coefficients_and_interactions, self.n_interaction_endmembers ) + self.n_W_mbr = len(endmember_coefficients_and_interactions) - self.transformation_matrix = transformation_matrix + if transformation_matrix is None: + self.dqdp = np.eye(self.n_endmembers) + else: + self.dqdp = transformation_matrix + self.n_transformed_endmembers = len(self.dqdp) self.reset() def reset(self): @@ -1091,6 +1172,7 @@ def set_composition(self, molar_fractions): """ self.reset() self.molar_fractions = molar_fractions + self.trans_fractions = np.einsum("ij, j->i", self.dqdp, molar_fractions) def set_state(self, pressure, temperature): """ @@ -1100,259 +1182,307 @@ def set_state(self, pressure, temperature): for mbr in self.interaction_endmembers: mbr.set_state(pressure, temperature) - def make_interaction_arrays(self, Ws): - n_Ws = len(Ws) - for i, W in enumerate(Ws): - if not all(W[i] <= W[i + 1] for i in range(3, len(W) - 1)): - raise Exception( - f"Interaction parameter {i+1}/{n_Ws} must be " - "upper triangular (i<=j<=k<=...<=z)" - ) - if not W[3] < W[-1]: - raise Exception( - f"Interaction parameter {i+1}/{n_Ws} must not lie on the " - "first diagonal (i=j=k=...=z)" - ) - - W_arrays = [] - - dims = sorted(list(set([len(W) - 3 for W in Ws]))) - for dim in dims: - coords = [ - [0, *W[3:]] for W in Ws if len(W) == dim + 3 and np.abs(W[0]) > 1.0e-10 - ] - coords.extend( - [ - [1, *W[3:]] - for W in Ws - if len(W) == dim + 3 and np.abs(W[1]) > 1.0e-10 - ] - ) - coords.extend( - [ - [2, *W[3:]] - for W in Ws - if len(W) == dim + 3 and np.abs(W[2]) > 1.0e-10 - ] + def _make_prefactors_and_exponents(self, a): + """ + This hidden method makes an Interaction object + containing convenience arrays to facilitate rapid + calculation of the first and second compositional + derivatives of the Gibbs, entropy and volume + excesses. + + :param a: A list of polynomial indices and exponents. + :type a: list + :return: An interaction object containing the + indices (inds), exponents (expts), + first derivative prefactors (j_r) and exponents (m_jr) and + second derivative prefactors (f_rs) and exponents (m_jrs). + :rtype: Interaction dataclass instance + """ + # assign exponent matrix for first derivatives + if len(a) % 2 != 0: + raise Exception( + "Interaction must have paired indices " f"and exponents, currently {a}" ) - coords = list(zip(*coords)) - - data = [W[0] for W in Ws if len(W) == dim + 3 and np.abs(W[0]) > 1.0e-10] - data.extend( - [W[1] for W in Ws if len(W) == dim + 3 and np.abs(W[1]) > 1.0e-10] + indices = a[::2] + exponents = np.array(a[1::2]) + n_indices = len(indices) + m_jr = np.zeros((n_indices, self.n_endmembers)) + m_jr[:, :] = exponents[:, np.newaxis] + m_jr[list(range(len(m_jr))), indices] -= np.ones(len(indices)) + + # assign factors for first derivatives + f_jr = np.zeros((n_indices, self.n_endmembers)) + f_jr[:, indices] = 1.0 + f_jr[range(len(m_jr)), indices] = exponents + f_r = np.prod(f_jr, axis=0) + + # assign exponent matrix for second derivatives + m_jrs = np.zeros((n_indices, self.n_endmembers, self.n_endmembers)) + m_jrs[:, :, :] = exponents[:, np.newaxis, np.newaxis] + + m_jrs -= ( + np.einsum( + "ij, k->ijk", np.eye(self.n_endmembers), np.ones(self.n_endmembers) ) - data.extend( - [W[2] for W in Ws if len(W) == dim + 3 and np.abs(W[2]) > 1.0e-10] + + np.einsum( + "ik, j->ijk", np.eye(self.n_endmembers), np.ones(self.n_endmembers) ) - - shape = [3] # First dimension is for E, S, V - shape.extend([self.n_endmembers for i in range(dim)]) - - shape = tuple(shape) - s = sparse.COO(coords, data, shape=shape).todense() - W_arrays.append((dim, s)) - return W_arrays - - def make_endmember_interaction_arrays(self, Ws): - n_Ws = len(Ws) - n_int = self.n_interaction_endmembers + )[indices] + + # assign factors for second derivatives + f_jrs = np.zeros((n_indices, self.n_endmembers, self.n_endmembers)) + f_jrs[np.ix_(range(n_indices), indices, indices)] = 1.0 + + ival = np.array([indices, exponents]).T + for k, (i, iexp) in enumerate(ival): + for l, (j, jexp) in enumerate(ival[k:]): + if i == j: + f_jrs[k, i, i] = iexp * (iexp - 1) + else: + f_jrs[k, i, j] = iexp + f_jrs[k, j, i] = iexp + f_jrs[l, j, i] = jexp + f_jrs[l, i, j] = jexp + + f_rs = np.prod(f_jrs, axis=0) + return Interaction(indices, exponents, f_r, m_jr, f_rs, m_jrs) + + def _make_interaction_arrays(self, Ws, n): + """ + A hidden convenience function that splits each + excess term into _summary_ + + :param Ws: List of interactions in form input by the user. + :type Ws: list of lists + :param n: Number of interaction values (e.g. 3 for ESV interactions) + :type n: integer + :return: The processed interaction values, indices and exponents. + :rtype: tuple of 2D numpy array and list of Interaction objects + """ + values = np.empty((len(Ws), n)) + interactions = [] for i, W in enumerate(Ws): - if not all(W[i] <= W[i + 1] for i in range(n_int, len(W) - 1)): - raise Exception( - f"Interaction parameter {i+1}/{n_Ws} must be " - f"upper triangular (i<=j<=k<=...<=z)\n(value = {W})" - ) - if not W[n_int] < W[-1]: - raise Exception( - f"Interaction parameter {i+1}/{n_Ws} must not lie on the " - f"first diagonal (i=j=k=...=z)\n(value = {W})" - ) - - W_arrays = [] - - dims = sorted(list(set([len(W) - n_int for W in Ws]))) - for dim in dims: - coords = [] - data = [] - for i in range(n_int): - coords.extend( - [ - [i, *W[n_int:]] - for W in Ws - if len(W) == dim + n_int and np.abs(W[i]) > 1.0e-10 - ] - ) - data.extend( - [ - W[i] - for W in Ws - if len(W) == dim + n_int and np.abs(W[i]) > 1.0e-10 - ] - ) - - coords = list(zip(*coords)) - - shape = [n_int] # First dimension is for the excess_endmembers - shape.extend([self.n_endmembers for i in range(dim)]) - - shape = tuple(shape) - s = sparse.COO(coords, data, shape=shape).todense() - W_arrays.append((dim, s)) - return W_arrays - - def transform_scalar_list(self, x, A): - if A is not None: - return A.dot(x) - else: - return x - - def transform_gradient_list(self, W, A): - if A is not None: - return np.einsum("ij, jk->ik", W, A) - else: - return W - - def transform_hessian_list(self, W, A): - if A is not None: - return np.einsum("ki, mij, lj->mkl", A, W, A) - else: - return W + values[i] = W[:n] + interactions.append(self._make_prefactors_and_exponents(W[n:])) + return values, interactions - def W_dots_x(self, dim, W_array, x): - if dim == 0: - return W_array - else: - ESVdim = [W_array] - ESVdim.extend([x for i in range(dim)]) - strng = string.ascii_lowercase[: W_array.ndim] - strng2 = ", ".join(strng[-dim:]) - strng += f", {strng2}" - return np.einsum(strng, *ESVdim) - - def rolled_array(self, W, istart): - W_array_rolled = deepcopy(W) - dim = len(W.shape) - 1 - - W_roll = deepcopy(W) - for i in range(dim - istart): - W_roll = np.moveaxis(W_roll, istart, -1) - W_array_rolled += W_roll - return W_array_rolled - - def get_scalar_list(self, x, W_arrays, A): - xp = self.transform_scalar_list(x, A) - xp_total = np.sum(xp) - ESV = np.zeros(W_arrays[0][-1].shape[0]) - for dim, W_array in W_arrays: - ESV += self.W_dots_x(dim, W_array, xp) / np.power(xp_total, dim - 1.0) - return ESV - - def get_gradient_list(self, x, W_arrays, A): - xp = self.transform_scalar_list(x, A) - xp_total = np.sum(xp) - - dESVdx = np.zeros(W_arrays[0][-1].shape[:2]) - - for dim, W_array in W_arrays: - dESVdx += ( - -(dim - 1) - / np.power(xp_total, dim) - * self.W_dots_x(dim, W_array, xp)[:, np.newaxis] + @material_property + def _c_xs(self): + """ + Calculates the product of molar fractions raised to the exponents in + each polynomial interaction. + :return: Interaction postfactors + (not including the values of the interactions themselves) + :rtype: 1D numpy array + """ + c = np.empty(self.n_W_ESV + self.n_W_mbr) + for i in range(self.n_W_ESV): + c[i] = np.prod( + np.power(self.trans_fractions[self.c_ESV[i].inds], self.c_ESV[i].expts), + axis=0, ) - - W_array_rolled = self.rolled_array(W_array, 1) - dESVdx += self.W_dots_x(dim - 1, W_array_rolled, xp) / np.power( - xp_total, dim - 1.0 + for i in range(self.n_W_mbr): + j = i + self.n_W_ESV + c[j] = np.prod( + np.power(self.trans_fractions[self.c_mbr[i].inds], self.c_mbr[i].expts), + axis=0, ) - return self.transform_gradient_list(dESVdx, A) - - def get_hessian_list(self, x, W_arrays, A): - xp = self.transform_scalar_list(x, A) - xp_total = np.sum(xp) - d2ESVdx2 = np.zeros(W_arrays[0][-1].shape[:3]) + return c - for dim, W_array in W_arrays: - d2ESVdx2 += ( - (dim - 1) - * dim - / np.power(xp_total, dim + 1) - * self.W_dots_x(dim, W_array, xp)[:, np.newaxis, np.newaxis] + @material_property + def _dc_xsdq(self): + """ + Calculates the first derivative of the product of + molar fractions raised to the exponents in each polynomial interaction + with respect to the transformed endmember proportions. + :return: Interaction postfactors + (not including the values of the interactions themselves) + :rtype: 2D numpy array + """ + c = np.empty((self.n_W_ESV + self.n_W_mbr, self.n_transformed_endmembers)) + for i in range(self.n_W_ESV): + c[i] = self.c_ESV[i].f_r * np.prod( + np.power( + self.trans_fractions[self.c_ESV[i].inds], self.c_ESV[i].m_jr.T + ), + axis=1, ) + for i in range(self.n_W_mbr): + j = i + self.n_W_ESV + c[j] = self.c_mbr[i].f_r * np.prod( + np.power( + self.trans_fractions[self.c_mbr[i].inds], self.c_mbr[i].m_jr.T + ), + axis=1, + ) + return c - W_array_rolled = self.rolled_array(W_array, 1) - - f = self.W_dots_x(dim - 1, W_array_rolled, xp) - h = f[:, np.newaxis, :] + f[:, :, np.newaxis] - - d2ESVdx2 += -(dim - 1) / np.power(xp_total, dim) * h + @material_property + def _d2c_xsdqdq(self): + """ + Calculates the second derivative of the product of + molar fractions raised to the exponents in each polynomial interaction + with respect to the transformed endmember proportions. + :return: Interaction postfactors + (not including the values of the interactions themselves) + :rtype: 3D numpy array + """ + c = np.empty( + ( + self.n_W_ESV + self.n_W_mbr, + self.n_transformed_endmembers, + self.n_transformed_endmembers, + ) + ) + for i in range(self.n_W_ESV): + c[i] = self.c_ESV[i].f_rs * np.prod( + np.power( + self.trans_fractions[self.c_ESV[i].inds], self.c_ESV[i].m_jrs.T + ), + axis=2, + ) + for i in range(self.n_W_mbr): + j = i + self.n_W_ESV + c[j] = self.c_mbr[i].f_rs * np.prod( + np.power( + self.trans_fractions[self.c_mbr[i].inds], self.c_mbr[i].m_jrs.T + ), + axis=2, + ) + return c - W_array_rolled_2 = self.rolled_array(W_array_rolled, 2) - g = self.W_dots_x(dim - 2, W_array_rolled_2, xp) + @material_property + def _dqdn(self): + """ + The first derivative of the transformed endmember proportions + with respect to the original endmember amounts. - d2ESVdx2 += g / np.power(xp_total, dim - 1.0) + :return: dqdn + :rtype: 2D numpy array + """ + return np.einsum( + "ik, km->im", + self.dqdp, + dpdn(self.molar_fractions, 1.0, self.ones, self.eye), + ) - return self.transform_hessian_list(d2ESVdx2, A) + @material_property + def _dc_xsdn(self): + """ + Calculates the first derivative of the product of + molar fractions raised to the exponents in each polynomial interaction + with respect to the endmember amounts. + :return: Interaction postfactors + (not including the values of the interactions themselves) + :rtype: 2D numpy array + """ + return np.einsum("ir, rm->im", self._dc_xsdq, self._dqdn) + np.einsum( + "i, m->im", self._c_xs, self.ones + ) @material_property - def ESV_scalar_list(self): - if self.W_ESV is None: - return np.zeros((3)) - else: - return self.get_scalar_list( - self.molar_fractions, self.W_ESV, self.transformation_matrix - ) + def _d2c_xsdndn(self): + """ + Calculates the second derivative of the product of + molar fractions raised to the exponents in each polynomial interaction + with respect to the endmember amounts. + :return: Interaction postfactors + (not including the values of the interactions themselves) + :rtype: 3D numpy array + """ + n = 1.0 + a1 = np.einsum("irs, sn, rm", self._d2c_xsdqdq, self._dqdn, self._dqdn) + a2a = 1.0 / n * (np.einsum("ir, rm, n", self._dc_xsdq, self._dqdn, self.ones)) + a2 = a2a + np.einsum("imn->inm", a2a) + a3 = np.einsum( + "ir, rk, knm->inm", + self._dc_xsdq, + self.dqdp, + d2pdndn(self.molar_fractions, n * n, self.ones, self.eyeones), + ) + return a1 + a2 + a3 @material_property - def ESV_gradient_list(self): + def _ESV_gradient_list(self): + """ + Calculates the first derivative of the + internal energy, entropy and volume + with respect to the endmember amounts. + + :return: dEdn, dSdn and dVdn. + :rtype: 2D numpy array + """ if self.W_ESV is None: return np.zeros((3, self.n_endmembers)) else: - return self.get_gradient_list( - self.molar_fractions, self.W_ESV, self.transformation_matrix - ) + return np.einsum("ij, ik->jk", self.W_ESV, self._dc_xsdn[: self.n_W_ESV]) @material_property - def ESV_hessian_list(self): + def _ESV_hessian_list(self): + """ + Calculates the second derivative of the + internal energy, entropy and volume + with respect to the endmember amounts. + + :return: d2Edndn, d2Sdndn and d2Vdndn. + :rtype: 3D numpy array + """ if self.W_ESV is None: return np.zeros((3, self.n_endmembers, self.n_endmembers)) else: - return self.get_hessian_list( - self.molar_fractions, self.W_ESV, self.transformation_matrix + return np.einsum( + "ij, ikl->jkl", self.W_ESV, self._d2c_xsdndn[: self.n_W_ESV] ) @material_property def mbr_scalar_list(self): + """ + :return: The value with which to multiply + each interaction endmember property to get the + total excess property. + :rtype: 1D numpy array + """ if self.W_mbr is None: return np.zeros((0)) else: - return self.get_scalar_list( - self.molar_fractions, self.W_mbr, self.transformation_matrix - ) + return np.einsum("ij, i->j", self.W_mbr, self._c_xs[self.n_W_ESV :]) @material_property - def mbr_gradient_list(self): + def _mbr_gradient_list(self): + """ + :return: The values with which to multiply + each interaction endmember property to get the + first derivatives of the total excess property + with respect to the endmember amounts. + :rtype: 2D numpy array + """ if self.W_mbr is None: return np.zeros((0, self.n_endmembers)) else: - return self.get_gradient_list( - self.molar_fractions, self.W_mbr, self.transformation_matrix - ) + return np.einsum("ij, ik->jk", self.W_mbr, self._dc_xsdn[self.n_W_ESV :]) @material_property - def mbr_hessian_list(self): + def _mbr_hessian_list(self): + """ + :return: The values with which to multiply + each interaction endmember property to get the + second derivatives of the total excess property + with respect to the endmember amounts. + :rtype: 3D numpy array + """ if self.W_mbr is None: return np.zeros((0, self.n_endmembers, self.n_endmembers)) else: - return self.get_hessian_list( - self.molar_fractions, self.W_mbr, self.transformation_matrix + return np.einsum( + "ij, ikl->jkl", self.W_mbr, self._d2c_xsdndn[self.n_W_ESV :] ) def _non_ideal_excess_partial_gibbs(self, pressure, temperature, molar_fractions): - dEdx, dSdx, dVdx = self.ESV_gradient_list - mbr_gradients = self.mbr_gradient_list + dEdn, dSdn, dVdn = self._ESV_gradient_list + mbr_gradients = self._mbr_gradient_list mbr_gibbs = np.array([mbr.gibbs for mbr in self.interaction_endmembers]) - gibbs = dEdx - temperature * dSdx + pressure * dVdx + gibbs = dEdn - temperature * dSdn + pressure * dVdn gibbs += np.einsum("i, ij->j", mbr_gibbs, mbr_gradients) return gibbs @@ -1372,55 +1502,55 @@ def excess_partial_entropies(self, pressure, temperature, molar_fractions): self, temperature, molar_fractions ) - dSdx = self.ESV_gradient_list[1] - mbr_gradients = self.mbr_gradient_list + dSdn = self._ESV_gradient_list[1] + mbr_gradients = self._mbr_gradient_list mbr_entropies = np.array([mbr.S for mbr in self.interaction_endmembers]) - non_ideal_entropies = dSdx + np.einsum("i, ij->j", mbr_entropies, mbr_gradients) + non_ideal_entropies = dSdn + np.einsum("i, ij->j", mbr_entropies, mbr_gradients) return ideal_entropies + non_ideal_entropies def excess_partial_volumes(self, pressure, temperature, molar_fractions): - dVdx = self.ESV_gradient_list[2] - mbr_gradients = self.mbr_gradient_list + dVdn = self._ESV_gradient_list[2] + mbr_gradients = self._mbr_gradient_list mbr_volumes = np.array([mbr.V for mbr in self.interaction_endmembers]) - return dVdx + np.einsum("i, ij->j", mbr_volumes, mbr_gradients) + return dVdn + np.einsum("i, ij->j", mbr_volumes, mbr_gradients) def gibbs_hessian(self, pressure, temperature, molar_fractions): ideal_entropy_hessian = IdealSolution._ideal_entropy_hessian( self, temperature, molar_fractions ) - d2Edx2, d2Sdx2, d2Vdx2 = self.ESV_hessian_list - mbr_hessian = self.mbr_hessian_list + d2Edndn, d2Sdndn, d2Vdndn = self._ESV_hessian_list + mbr_hessian = self._mbr_hessian_list mbr_gibbs = np.array([mbr.gibbs for mbr in self.interaction_endmembers]) - d2Gdx2 = d2Edx2 - temperature * d2Sdx2 + pressure * d2Vdx2 - d2Gdx2 += np.einsum("i, ijk->jk", mbr_gibbs, mbr_hessian) + d2Gdndn = d2Edndn - temperature * d2Sdndn + pressure * d2Vdndn + d2Gdndn += np.einsum("i, ijk->jk", mbr_gibbs, mbr_hessian) - return d2Gdx2 - temperature * ideal_entropy_hessian + return d2Gdndn - temperature * ideal_entropy_hessian def entropy_hessian(self, pressure, temperature, molar_fractions): ideal_entropy_hessian = IdealSolution._ideal_entropy_hessian( self, temperature, molar_fractions ) - d2Sdx2 = self.ESV_hessian_list[1] - mbr_hessian = self.mbr_hessian_list + d2Sdndn = self._ESV_hessian_list[1] + mbr_hessian = self._mbr_hessian_list mbr_entropies = np.array([mbr.S for mbr in self.interaction_endmembers]) - d2Sdx2 += np.einsum("i, ijk->jk", mbr_entropies, mbr_hessian) + d2Sdndn += np.einsum("i, ijk->jk", mbr_entropies, mbr_hessian) - return d2Sdx2 + ideal_entropy_hessian + return d2Sdndn + ideal_entropy_hessian def volume_hessian(self, pressure, temperature, molar_fractions): - d2Vdx2 = self.ESV_hessian_list[2] - mbr_hessian = self.mbr_hessian_list + d2Vdndn = self._ESV_hessian_list[2] + mbr_hessian = self._mbr_hessian_list mbr_volumes = np.array([mbr.V for mbr in self.interaction_endmembers]) - d2Vdx2 += np.einsum("i, ijk->jk", mbr_volumes, mbr_hessian) + d2Vdndn += np.einsum("i, ijk->jk", mbr_volumes, mbr_hessian) - return d2Vdx2 + return d2Vdndn def Cp_excess(self): mbr_scalar = self.mbr_scalar_list diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py index c3a48643..a1354714 100644 --- a/tests/test_solidsolution.py +++ b/tests/test_solidsolution.py @@ -20,6 +20,7 @@ from burnman.minerals.SLB_2011 import mg_post_perovskite from burnman.minerals.SLB_2011 import fe_post_perovskite from burnman.minerals.SLB_2011 import al_post_perovskite +from burnman.tools.eos import check_eos_consistency class forsterite(Mineral): @@ -247,23 +248,23 @@ def __init__(self, molar_fractions=None): ], ESV_interactions=[ # Wj and Wi terms in Helffrich and Wood - [10.0e3, 1.0, 0.0, 0, 1, 1], - [-10.0e3, -2.0, 0.0, 0, 0, 1], - [5.0e3, 0.0, 0.0, 0, 2, 2], - [3.0e3, 1.0, 0.0, 0, 0, 2], - [-10.0e3, 0.0, 0.0, 1, 2, 2], - [-10.0e3, 0.0, 0.0, 1, 1, 2], + [10.0e3, 1.0, 0.0, 0, 1, 1, 2], + [-10.0e3, -2.0, 0.0, 0, 2, 1, 1], + [5.0e3, 0.0, 0.0, 0, 1, 2, 2], + [3.0e3, 1.0, 0.0, 0, 2, 2, 1], + [-10.0e3, 0.0, 0.0, 1, 1, 2, 2], + [-10.0e3, 0.0, 0.0, 1, 2, 2, 1], # 0.5*Wk terms in Helffrich and Wood # Because there are only three endmembers, there # is only one additional term for each term above: - [5.0e3, 0.5, 0.0, 0, 1, 2], - [-5.0e3, -1.0, 0.0, 0, 1, 2], - [2.5e3, 0.0, 0.0, 0, 1, 2], - [1.5e3, 0.5, 0.0, 0, 1, 2], - [-5.0e3, 0.0, 0.0, 0, 1, 2], - [-5.0e3, 0.0, 0.0, 0, 1, 2], + [5.0e3, 0.5, 0.0, 0, 1, 1, 1, 2, 1], + [-5.0e3, -1.0, 0.0, 0, 1, 1, 1, 2, 1], + [2.5e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], + [1.5e3, 0.5, 0.0, 0, 1, 1, 1, 2, 1], + [-5.0e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], + [-5.0e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], # Ternary term - [3.0e3, 0.0, 0.0, 0, 1, 2], + [3.0e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], ], ) @@ -288,23 +289,46 @@ def __init__(self, molar_fractions=None): ], ESV_interactions=[ # Wj and Wi terms in Helffrich and Wood - [10.0e3, 1.0, 0.0, 0, 1, 1], - [-10.0e3, -2.0, 0.0, 0, 0, 1], - [5.0e3, 0.0, 0.0, 0, 2, 2], - [3.0e3, 1.0, 0.0, 0, 0, 2], - [-10.0e3, 0.0, 0.0, 1, 2, 2], - [-10.0e3, 0.0, 0.0, 1, 1, 2], + [10.0e3, 1.0, 0.0, 0, 1, 1, 2], + [-10.0e3, -2.0, 0.0, 0, 2, 1, 1], + [5.0e3, 0.0, 0.0, 0, 1, 2, 2], + [3.0e3, 1.0, 0.0, 0, 2, 2, 1], + [-10.0e3, 0.0, 0.0, 1, 1, 2, 2], + [-10.0e3, 0.0, 0.0, 1, 2, 2, 1], # 0.5*Wk terms in Helffrich and Wood # Because there are only three endmembers, there # is only one additional term for each term above: - [5.0e3, 0.5, 0.0, 0, 1, 2], - [-5.0e3, -1.0, 0.0, 0, 1, 2], - [2.5e3, 0.0, 0.0, 0, 1, 2], - [1.5e3, 0.5, 0.0, 0, 1, 2], - [-5.0e3, 0.0, 0.0, 0, 1, 2], - [-5.0e3, 0.0, 0.0, 0, 1, 2], + [5.0e3, 0.5, 0.0, 0, 1, 1, 1, 2, 1], + [-5.0e3, -1.0, 0.0, 0, 1, 1, 1, 2, 1], + [2.5e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], + [1.5e3, 0.5, 0.0, 0, 1, 1, 1, 2, 1], + [-5.0e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], + [-5.0e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], # Ternary term - [3.0e3, 0.0, 0.0, 0, 1, 2], + [3.0e3, 0.0, 0.0, 0, 1, 1, 1, 2, 1], + ], + transformation_matrix=np.array( + [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]] + ), + ) + + burnman.SolidSolution.__init__(self, molar_fractions) + + +class two_site_ss_polynomial_high_order_transformed(burnman.SolidSolution): + def __init__(self, molar_fractions=None): + self.name = "two_site_ss (high order from polynomial)" + self.solution_model = PolynomialSolution( + endmembers=[ + [forsterite(), "[Mg]3[Mg1/2Si1/2]2Si3O12"], + [forsterite(), "[Mg]3[Al]2Si3O12"], + [forsterite(), "[Fe]3[Al]2Si3O12"], + ], + ESV_interactions=[ + # Wj and Wi terms in Helffrich and Wood + [10.0e3, 1.0, 0.0, 0, 1, 1, 2], + [-10.0e3, -2.0, 0.0, 0, 2, 1, 6], + [5.0e3, 0.5, 0.0, 0, 6, 1, 1, 2, 1], ], transformation_matrix=np.array( [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]] @@ -325,14 +349,30 @@ def __init__(self, molar_fractions=None): class polynomial_mbr_mixing(burnman.SolidSolution): def __init__(self, molar_fractions=None): - self.name = "binary mbr mixing" + self.name = "binary mbr mixing (polynomial)" self.solution_model = PolynomialSolution( endmembers=[ [forsterite(), "[Mg]2"], [fayalite(), "[Mg]2"], ], interaction_endmembers=[fo2, fo1], - endmember_coefficients_and_interactions=[[1.0, -1.0, 0, 1]], + endmember_coefficients_and_interactions=[[1.0, -1.0, 0, 1, 1, 1]], + ) + + burnman.SolidSolution.__init__(self, molar_fractions) + + +class regular_copy_polynomial_mbr_mixing(burnman.SolidSolution): + def __init__(self, molar_fractions=None): + self.name = "binary mbr mixing (regular)" + self.solution_model = SymmetricRegularSolution( + endmembers=[ + [forsterite(), "[Mg]2"], + [fayalite(), "[Mg]2"], + ], + energy_interaction=[[1200.0]], + entropy_interaction=[[5.0]], + volume_interaction=[[1.0e-6]], ) burnman.SolidSolution.__init__(self, molar_fractions) @@ -352,9 +392,25 @@ def __init__(self, molar_fractions=None): [forsterite(), "[Mg]2"], [fayalite(), "[Mg]2"], ], - ESV_interactions=[[600.0, 2.5, 0.5e-6, 0, 1]], + ESV_interactions=[[600.0, 2.5, 0.5e-6, 0, 1, 1, 1]], interaction_endmembers=[fo3, fo1], - endmember_coefficients_and_interactions=[[1.0, -1.0, 0, 1]], + endmember_coefficients_and_interactions=[[1.0, -1.0, 0, 1, 1, 1]], + ) + + burnman.SolidSolution.__init__(self, molar_fractions) + + +class polynomial_esv_and_mbr_mixing_2(burnman.SolidSolution): + def __init__(self, molar_fractions=None): + self.name = "two_site_ss (subregular symmetric)" + self.solution_model = PolynomialSolution( + endmembers=[ + [forsterite(), "[Mg]2"], + [fayalite(), "[Mg]2"], + ], + ESV_interactions=[[600.0, 2.5, 0.5e-6, 0, 1, 1, 1]], + interaction_endmembers=[fayalite(), fo1], + endmember_coefficients_and_interactions=[[1.0, -1.0, 0, 1, 1, 1]], ) burnman.SolidSolution.__init__(self, molar_fractions) @@ -754,6 +810,61 @@ def test_subregular_model_ternary_hessian_multicomponent_change(self): self.assertArraysAlmostEqual(H0.dot(df), dGdx2 - dGdx1) + def test_polynomial_model_partial_entropy_multicomponent_change(self): + ss = two_site_ss_polynomial_ternary() + f0 = np.array([0.25, 0.35, 0.4]) + ss.set_composition(f0) + ss.set_state(1.0e9, 1000.0) + + dSdx = ss.partial_entropies + + df = 0.0001 + dSdx2 = np.empty(3) + for i, f_mod in enumerate(np.eye(3) * df): + ss.set_composition((f0 - f_mod / 2) / (1.0 - df / 2.0)) + S0 = ss.S * (1.0 - df / 2.0) + ss.set_composition((f0 + f_mod / 2) / (1.0 + df / 2.0)) + S1 = ss.S * (1.0 + df / 2.0) + + dSdx2[i] = (S1 - S0) / df + + self.assertArraysAlmostEqual(dSdx, dSdx2) + + def test_polynomial_high_order_transformed(self): + ss = two_site_ss_polynomial_high_order_transformed() + f0 = [0.25, 0.35, 0.4] + ss.set_composition(f0) + ss.set_state(1.0e5, 300.0) + HG0 = ss.gibbs_hessian + HS0 = ss.entropy_hessian + HV0 = ss.volume_hessian + dGdx0 = ss.partial_gibbs + dSdx0 = ss.partial_entropies + dVdx0 = ss.partial_volumes + + df = np.array([2.0e-6, -1.5e-6, -0.5e-6]) + ss.set_composition(f0 - df / 2.0) + G1 = ss.gibbs + S1 = ss.S + V1 = ss.V + dGdx1 = ss.partial_gibbs + dSdx1 = ss.partial_entropies + dVdx1 = ss.partial_volumes + ss.set_composition(f0 + df / 2.0) + G2 = ss.gibbs + S2 = ss.S + V2 = ss.V + dGdx2 = ss.partial_gibbs + dSdx2 = ss.partial_entropies + dVdx2 = ss.partial_volumes + + self.assertArraysAlmostEqual( + [dGdx0.dot(df), dSdx0.dot(df), dVdx0.dot(df)], [G2 - G1, S2 - S1, V2 - V1] + ) + self.assertArraysAlmostEqual(HG0.dot(df), dGdx2 - dGdx1) + self.assertArraysAlmostEqual(HS0.dot(df), dSdx2 - dSdx1) + self.assertArraysAlmostEqual(HV0.dot(df), dVdx2 - dVdx1) + def test_polynomial_using_subregular_scalars(self): ss1 = two_site_ss_subregular_ternary() ss2 = two_site_ss_polynomial_ternary() @@ -830,18 +941,31 @@ def test_polynomial_using_subregular_volume(self): self.assertArraysAlmostEqual( ss1.volume_hessian.flatten(), ss2.volume_hessian.flatten() ) + self.assertArraysAlmostEqual(ss1.partial_gibbs, ss2.partial_gibbs) + self.assertArraysAlmostEqual(ss1.partial_entropies, ss2.partial_entropies) + self.assertArraysAlmostEqual(ss1.partial_volumes, ss2.partial_volumes) + + def test_polynomial_consistency(self): + ss = polynomial_esv_and_mbr_mixing_2() + ss.set_composition([0.4, 0.6]) + self.assertTrue(check_eos_consistency(ss)) def test_polynomial_mbr_mixing(self): - ss = polynomial_mbr_mixing() + ss1 = polynomial_mbr_mixing() + ss2 = regular_copy_polynomial_mbr_mixing() f0 = [0.25, 0.75] - ss.set_state(1.0e5, 300.0) - ss.set_composition(f0) + for ss in [ss1, ss2]: + ss.set_state(1.0e5, 300.0) + ss.set_composition(f0) self.assertAlmostEqual( - ss.excess_gibbs, 0.25 * 0.75 * (1200.0 - 300.0 * 5 + 1.0e5 * 1.0e-6) + ss1.excess_gibbs, 0.25 * 0.75 * (1200.0 - 300.0 * 5 + 1.0e5 * 1.0e-6) ) - self.assertAlmostEqual(ss.excess_entropy, 0.25 * 0.75 * 5) - self.assertAlmostEqual(ss.excess_volume, 0.25 * 0.75 * 1.0e-6) + self.assertAlmostEqual(ss1.excess_entropy, 0.25 * 0.75 * 5) + self.assertAlmostEqual(ss1.excess_volume, 0.25 * 0.75 * 1.0e-6) + self.assertArraysAlmostEqual(ss1.partial_gibbs, ss2.partial_gibbs) + self.assertArraysAlmostEqual(ss1.partial_entropies, ss2.partial_entropies) + self.assertArraysAlmostEqual(ss1.partial_volumes, ss2.partial_volumes) def test_polynomial_ESV_and_mbr_mixing(self): ss = polynomial_esv_and_mbr_mixing() @@ -860,10 +984,19 @@ def test_polynomial_ternary_transformed(self): for ss in [ss1, ss2]: ss.set_state(1.0e5, 300.0) - ss1.set_composition([0.25, 0.35, 0.4]) - ss2.set_composition([0.4, 0.25, 0.35]) - + f = np.array([0.25, 0.35, 0.4]) + ss1.set_composition(f) + # roll forward to set composition + ss2.set_composition(np.roll(f, 1)) self.assertAlmostEqual(ss1.excess_gibbs, ss2.excess_gibbs) + # roll back to compare properties + self.assertArraysAlmostEqual(ss1.partial_gibbs, np.roll(ss2.partial_gibbs, -1)) + self.assertArraysAlmostEqual( + ss1.partial_entropies, np.roll(ss2.partial_entropies, -1) + ) + self.assertArraysAlmostEqual( + ss1.partial_volumes, np.roll(ss2.partial_volumes, -1) + ) def test_subregular(self): ss0 = two_site_ss()