diff --git a/.github/workflows/pytest.yaml b/.github/workflows/pytest.yaml index b8c705f..c291ad0 100644 --- a/.github/workflows/pytest.yaml +++ b/.github/workflows/pytest.yaml @@ -11,7 +11,7 @@ jobs: fail-fast: false max-parallel: 100 matrix: - python-version: [3.7, 3.8, 3.9] + python-version: [3.8, 3.9] pycalphad_develop_version: [true, false] steps: diff --git a/scheil/simulate.py b/scheil/simulate.py index beba5d1..3202968 100644 --- a/scheil/simulate.py +++ b/scheil/simulate.py @@ -1,5 +1,6 @@ import sys import numpy as np +from collections import defaultdict from pycalphad import equilibrium, variables as v from pycalphad.codegen.callables import build_phase_records from pycalphad.core.calculate import _sample_phase_constitution @@ -107,10 +108,17 @@ def simulate_scheil_solidification(dbf, comps, phases, composition, solid_phases = sorted((set(phases) | filtered_disordered_phases) - {liquid_phase_name}) temp = start_temperature independent_comps = sorted([str(comp)[2:] for comp in composition.keys()]) + x_liquid = {comp: [composition[v.X(comp)]] for comp in independent_comps} fraction_solid = [0.0] temperatures = [temp] phase_amounts = {ph: [0.0] for ph in solid_phases} + x_phases = defaultdict(dict) + Y_phases = defaultdict(dict) + for ph in solid_phases: + Y_phases[ph] = [np.nan] + for comp in independent_comps: + x_phases[ph][comp] = [np.nan] if adaptive: dof_dict = {phase_name: list(map(len, mod.constituents)) for phase_name, mod in models.items()} @@ -183,12 +191,22 @@ def simulate_scheil_solidification(dbf, comps, phases, composition, for solid_phase in solid_phases: if solid_phase not in eq_phases: phase_amounts[solid_phase].append(0.0) + Y_phases[solid_phase].append(np.nan) + for comp in independent_comps: + x_phases[solid_phase][comp].append(np.nan) continue np_tieline = np.nansum(eq.isel(vertex=eq_phases.index(solid_phase))["NP"].values.squeeze()) found_phase_amounts.append((solid_phase, np_tieline)) delta_fraction_solid = (1 - current_fraction_solid) * np_tieline current_fraction_solid += delta_fraction_solid phase_amounts[solid_phase].append(delta_fraction_solid) + vertex = sorted(np.nonzero(eq["Phase"].values.squeeze().flat == solid_phase))[0] + Y_phases[solid_phase].append(list(eq["Y"].isel(vertex=vertex).squeeze().values)) + for comp in independent_comps: + if np.sum(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)==0: + x_phases[solid_phase][comp].append(np.nan) + else: + x_phases[solid_phase][comp].append(float(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)) fraction_solid.append(current_fraction_solid) temperatures.append(temp) NL = 1 - fraction_solid[-1] @@ -217,11 +235,22 @@ def simulate_scheil_solidification(dbf, comps, phases, composition, for solid_phase in solid_phases: if solid_phase in eq_phases: amount = np.nansum(eq.isel(vertex=eq_phases.index(solid_phase))["NP"].values.squeeze()) - phase_amounts[solid_phase].append(float(amount) * (1 - current_fraction_solid)) + try: + phase_amounts[solid_phase].append(float(amount) * (1 - current_fraction_solid)) + except UnboundLocalError: + phases_solid = set(eq_phases) - {''} + raise ValueError("The input temperature is too low to support a liquid phase, it only contains "+', '.join(phases_solid)) else: phase_amounts[solid_phase].append(0.0) + vertex = sorted(np.nonzero(eq["Phase"].values.squeeze().flat == solid_phase))[0] + Y_phases[solid_phase].append(list(eq["Y"].isel(vertex=vertex).squeeze().values)) + for comp in independent_comps: + if np.sum(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)==0: + x_phases[solid_phase][comp].append(np.nan) + else: + x_phases[solid_phase][comp].append(float(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)) - return SolidificationResult(x_liquid, fraction_solid, temperatures, phase_amounts, converged, "scheil") + return SolidificationResult(x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, "scheil") def simulate_equilibrium_solidification(dbf, comps, phases, composition, @@ -291,6 +320,12 @@ def simulate_equilibrium_solidification(dbf, comps, phases, composition, fraction_solid = [] phase_amounts = {ph: [] for ph in solid_phases} # instantaneous phase amounts cum_phase_amounts = {ph: [] for ph in solid_phases} + x_phases = defaultdict(dict) + Y_phases = defaultdict(dict) + for ph in solid_phases: + Y_phases[ph] = [] + for comp in independent_comps: + x_phases[ph][comp] = [] converged = False current_T = start_temperature if verbose: @@ -368,8 +403,16 @@ def simulate_equilibrium_solidification(dbf, comps, phases, composition, phase_amounts[solid_phase].append(amount) else: phase_amounts[solid_phase].append(amount - cum_phase_amounts[solid_phase][-2]) + vertex = sorted(np.nonzero(eq["Phase"].values.squeeze().flat == solid_phase))[0] + Y_phases[solid_phase].append(list(eq["Y"].isel(vertex=vertex).squeeze().values)) + for comp in independent_comps: + if np.sum(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)==0: + x_phases[solid_phase][comp].append(np.nan) + else: + x_phases[solid_phase][comp].append(float(eq["X"].isel(vertex=vertex).squeeze().sel(component=comp).values)) current_fraction_solid += amount + fraction_solid.append(current_fraction_solid) converged = True if np.isclose(fraction_solid[-1], 1.0) else False - return SolidificationResult(x_liquid, fraction_solid, temperatures, phase_amounts, converged, "equilibrium") + return SolidificationResult(x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, "equilibrium") diff --git a/scheil/solidification_result.py b/scheil/solidification_result.py index a070f67..25a52d6 100644 --- a/scheil/solidification_result.py +++ b/scheil/solidification_result.py @@ -38,8 +38,10 @@ class SolidificationResult(): """ - def __init__(self, x_liquid, fraction_solid, temperatures, phase_amounts, converged, method): + def __init__(self, x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, method): self.x_liquid = x_liquid + self.x_phases = x_phases + self.Y_phases = Y_phases self.fraction_solid = fraction_solid self.fraction_liquid = (1.0 - np.array(fraction_solid)).tolist() self.temperatures = temperatures @@ -57,6 +59,8 @@ def __repr__(self): def to_dict(self): d = { 'x_liquid': self.x_liquid, + 'x_phases': self.x_phases, + 'Y_phases': self.Y_phases, 'fraction_solid': self.fraction_solid, 'temperatures': self.temperatures, 'phase_amounts': self.phase_amounts, @@ -68,9 +72,11 @@ def to_dict(self): @classmethod def from_dict(cls, d): x_liquid = d['x_liquid'] + x_phases = d['x_phases'] + Y_phases = d['Y_phases'] fraction_solid = d['fraction_solid'] temperatures = d['temperatures'] phase_amounts = d['phase_amounts'] converged = d['converged'] method = d['method'] - return cls(x_liquid, fraction_solid, temperatures, phase_amounts, converged, method) + return cls(x_liquid, x_phases, Y_phases, fraction_solid, temperatures, phase_amounts, converged, method) diff --git a/tests/test_scheil_solidification.py b/tests/test_scheil_solidification.py index 8b2eb95..cd13fcb 100644 --- a/tests/test_scheil_solidification.py +++ b/tests/test_scheil_solidification.py @@ -26,6 +26,9 @@ def test_scheil_solidification_result_properties(): assert num_temperatures == len(sol_res.fraction_solid) assert all([num_temperatures == len(nphase) for nphase in sol_res.phase_amounts.values()]) assert all([num_temperatures == len(liq_comps) for liq_comps in sol_res.x_liquid.values()]) + for phase in sol_res.x_phases: + assert all([num_temperatures == len(phase_comps) for phase_comps in sol_res.x_phases[phase].values()]) + assert all([num_temperatures == len(phase_sf) for phase_sf in sol_res.Y_phases.values()]) assert all([num_temperatures == len(nphase) for nphase in sol_res.cum_phase_amounts.values()]) # final phase amounts are correct @@ -48,6 +51,8 @@ def test_scheil_solidification_result_properties(): assert rnd_trip_sol_res.fraction_liquid == sol_res.fraction_liquid assert rnd_trip_sol_res.fraction_solid == sol_res.fraction_solid assert rnd_trip_sol_res.x_liquid == sol_res.x_liquid + assert rnd_trip_sol_res.x_phases == sol_res.x_phases + assert rnd_trip_sol_res.Y_phases == sol_res.Y_phases assert rnd_trip_sol_res.cum_phase_amounts == sol_res.cum_phase_amounts assert rnd_trip_sol_res.phase_amounts == sol_res.phase_amounts assert rnd_trip_sol_res.temperatures == sol_res.temperatures