diff --git a/burnman/classes/material.py b/burnman/classes/material.py index 19558677..fd32de37 100644 --- a/burnman/classes/material.py +++ b/burnman/classes/material.py @@ -259,7 +259,7 @@ def unroll(self): """ raise NotImplementedError("need to implement unroll() in derived class!") - def evaluate(self, vars_list, pressures, temperatures): + def evaluate(self, vars_list, pressures, temperatures, molar_fractions=None): """ Returns an array of material properties requested through a list of strings at given pressure and temperature conditions. @@ -275,9 +275,11 @@ def evaluate(self, vars_list, pressures, temperatures): :param temperatures: ndlist or ndarray of float of temperatures in [K]. :type temperatures: :class:`numpy.array`, n-dimensional - :returns: Array returning all variables at given pressure/temperature values. + :returns: List or array returning all variables at given pressure/temperature values. output[i][j] is property vars_list[j] for temperatures[i] and pressures[i]. - :rtype: :class:`numpy.array`, n-dimensional + Attempts to return an array, falls back to a list if the returned properties + have different shapes. + :rtype: list or :class:`numpy.array`, n-dimensional """ old_pressure = self.pressure old_temperature = self.temperature @@ -286,11 +288,93 @@ def evaluate(self, vars_list, pressures, temperatures): assert pressures.shape == temperatures.shape - output = np.empty((len(vars_list),) + pressures.shape) + if molar_fractions is not None: + molar_fractions = np.array(molar_fractions) + assert temperatures.shape == molar_fractions.shape[:-1] + + # First, check the output types of all the requested variables: + self.set_state(pressures.flat[0], temperatures.flat[0]) + + output = [] + for j in range(len(vars_list)): + try: + var_shape = getattr(self, vars_list[j]).shape + except AttributeError: + var_shape = () + output.append(np.empty(pressures.shape + var_shape)) + for i, p in np.ndenumerate(pressures): + if molar_fractions is not None: + self.set_composition(molar_fractions[i]) self.set_state(p, temperatures[i]) for j in range(len(vars_list)): - output[(j,) + i] = getattr(self, vars_list[j]) + output[j][i] = getattr(self, vars_list[j]) + if old_pressure is None or old_temperature is None: + # do not set_state if old values were None. Just reset to None + # manually + self._pressure = self._temperature = None + self.reset() + else: + self.set_state(old_pressure, old_temperature) + + try: + output = np.array(output) + except ValueError: # if the lists are different shapes + pass + return output + + def evaluate_with_volumes( + self, vars_list, volumes, temperatures, molar_fractions=None + ): + """ + Returns an array of material properties requested through a list of strings + at given volume and temperature conditions. + At the end it resets the set_state to the original values. + The user needs to call set_method() before. + + :param vars_list: Variables to be returned for given conditions + :type vars_list: list of strings + + :param volumes: ndlist or ndarray of float of volumes in [m^3]. + :type volumes: :class:`numpy.array`, n-dimensional + + :param temperatures: ndlist or ndarray of float of temperatures in [K]. + :type temperatures: :class:`numpy.array`, n-dimensional + + :returns: List or array returning all variables at given pressure/temperature values. + output[i][j] is property vars_list[j] for temperatures[i] and pressures[i]. + Attempts to return an array, falls back to a list if the returned properties + have different shapes. + :rtype: list or :class:`numpy.array`, n-dimensional + """ + old_pressure = self.pressure + old_temperature = self.temperature + volumes = np.array(volumes) + temperatures = np.array(temperatures) + + assert volumes.shape == temperatures.shape + + if molar_fractions is not None: + molar_fractions = np.array(molar_fractions) + assert temperatures.shape == molar_fractions.shape[:-1] + + # First, check the output types of all the requested variables: + self.set_state_with_volume(volumes.flat[0], temperatures.flat[0]) + + output = [] + for j in range(len(vars_list)): + try: + var_shape = getattr(self, vars_list[j]).shape + except AttributeError: + var_shape = () + output.append(np.empty(volumes.shape + var_shape)) + + for i, v in np.ndenumerate(volumes): + if molar_fractions is not None: + self.set_composition(molar_fractions[i]) + self.set_state_with_volume(v, temperatures[i]) + for j in range(len(vars_list)): + output[j][i] = getattr(self, vars_list[j]) if old_pressure is None or old_temperature is None: # do not set_state if old values were None. Just reset to None # manually @@ -299,6 +383,10 @@ def evaluate(self, vars_list, pressures, temperatures): else: self.set_state(old_pressure, old_temperature) + try: + output = np.array(output) + except ValueError: # if the lists are different shapes + pass return output @property diff --git a/tests/test_averaging.py b/tests/test_averaging.py index a257ccd3..7cd3bcde 100644 --- a/tests/test_averaging.py +++ b/tests/test_averaging.py @@ -175,10 +175,10 @@ def test_non_present_non_rigid_phase(self): with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") rho, v_p, v_s, v_phi, K, G = rock.evaluate( - ["rho", "v_p", "v_s", "v_phi", "K_S", "G"], [10.0e9], [300.0] + ["rho", "v_p", "v_s", "v_phi", "K_S", "G"], 10.0e9, 300.0 ) - assert len(w) == 1 # we expect one error to be thrown when p_nonrigid = 0 - self.assertFloatEqual(150.901, G[0] / 1.0e9) + assert len(w) == 2 # we expect two errors to be thrown when p_nonrigid = 0 + self.assertFloatEqual(150.901, G / 1.0e9) def test_present_non_rigid_phase(self): rock = burnman.Composite([mypericlase(), my_nonrigid_mineral()], [0.5, 0.5]) @@ -186,10 +186,12 @@ def test_present_non_rigid_phase(self): with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always") rho, v_p, v_s, v_phi, K, G = rock.evaluate( - ["rho", "v_p", "v_s", "v_phi", "K_S", "G"], [10.0e9], [300.0] + ["rho", "v_p", "v_s", "v_phi", "K_S", "G"], 10.0e9, 300.0 ) - assert len(w) == 2 # we expect two errors to be thrown when p_nonrigid != 0 - self.assertFloatEqual(0.0, G[0] / 1.0e9) + assert ( + len(w) == 4 + ) # we expect four errors to be thrown when p_nonrigid != 0 + self.assertFloatEqual(0.0, G / 1.0e9) class Voigt(BurnManTest): diff --git a/tests/test_solidsolution.py b/tests/test_solidsolution.py index a1354714..73683975 100644 --- a/tests/test_solidsolution.py +++ b/tests/test_solidsolution.py @@ -1203,6 +1203,21 @@ def test_function_solution_low_proportions(self): ss[0].activity_coefficients, ss[1].activity_coefficients ) + def test_evaluate(self): + ol = olivine_ideal_ss() + molar_fractions = np.array([[0.1, 0.9], [0.2, 0.8]]) + pressures = [1.0e9, 2.0e9] + temperatures = [300.0, 600.0] + G1 = ol.evaluate(["gibbs"], pressures, temperatures, molar_fractions)[0] + + G2 = [] + for i in range(2): + ol.set_state(pressures[i], temperatures[i]) + ol.set_composition(molar_fractions[i]) + G2.append(ol.gibbs) + + self.assertArraysAlmostEqual(G1, G2) + if __name__ == "__main__": unittest.main()