diff --git a/burnman/tools/polytope.py b/burnman/tools/polytope.py index b4b833e8..c8471679 100644 --- a/burnman/tools/polytope.py +++ b/burnman/tools/polytope.py @@ -179,6 +179,10 @@ def simplify_composite_with_composition(composite, composition): solutions, this function will return a composite that only contains the Mg-bearing endmembers. + If the solutions have endmember proportions that are consistent with the + bulk composition, the site occupancies of the new solution models are set to + the values in the original models. + :param composite: The initial Composite object. :type composite: :class:`burnman.Composite` object @@ -239,7 +243,24 @@ def simplify_composite_with_composition(composite, composition): ) composite_changed = True - soln = transform_solution_to_new_basis(ph, new_basis) + + # If the composition is set and + # consistent with the new basis, + # make a new solution with the composition + # already set. + try: + sol = np.linalg.lstsq( + new_basis.T, ph.molar_fractions, rcond=None + ) + if sol[1][0] > 1.0e-10: + f = None + else: + f = sol[0] + except AttributeError: + f = None + soln = transform_solution_to_new_basis( + ph, new_basis, molar_fractions=f + ) new_phases.append(soln) else: logging.info( diff --git a/tests/test_polytope.py b/tests/test_polytope.py index ec52fbfe..d6b89d50 100644 --- a/tests/test_polytope.py +++ b/tests/test_polytope.py @@ -47,6 +47,21 @@ def test_simplify_composite(self): self.assertEqual(strings[0], "[Fe]3[Mg][Si]") self.assertEqual(strings[1], "[Mg]3[Mg][Si]") + def test_simplify_composite_and_composition(self): + gt = SLB_2011.garnet() + gt.set_composition([-0.1, 0.1, 0.0, 1.0, 0.0]) + ol = SLB_2011.mg_fe_olivine() + assemblage = Composite([ol, gt], [0.7, 0.3]) + composition = {"Fe": 3.0, "Mg": 1.0, "Si": 3.9, "O": 11.8} + + new_assemblage = simplify_composite_with_composition(assemblage, composition) + new_gt = new_assemblage.phases[1] + self.assertTrue(new_gt.n_endmembers == 2) + strings = list(sorted([e[1] for e in new_gt.endmembers])) + self.assertEqual(strings[0], "[Fe]3[Mg][Si]") + self.assertEqual(strings[1], "[Mg]3[Mg][Si]") + self.assertArraysAlmostEqual([0.1, 0.9], new_gt.molar_fractions) + if __name__ == "__main__": unittest.main()