preserve composition of solutions when simplifying composites

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(

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