Merge pull request #581 from bobmyhill/evaluate_composition

allow evaluation with varying composition

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

tests/test_averaging.py

@@ -175,21 +175,23 @@ 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])
         rock.set_averaging_scheme(avg.Reuss())
         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):

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()