Holding Point Push for sync

Author: amr8004

SStoTi64.ipynb

@@ -0,0 +1,111 @@
+"""This module is a script that reads the AMMap configuration YAML file and if it detects ``type: LC density`` under 
+the ``constraints`` key, it will proceed further. It will read the ``name`` key, append it with truncated (6 
+characters) hash of the YAML file, append it with list of elements forming the elemental spcae, and use it as the 
+name and if the directory with the name does not exist under the ``ammap/callables`` directory, it will create one 
+and also create a ``__init__.py`` file to make it a package. It will then read the ``ammap/templates/LCdensity.py``, 
+prepend it with constants based on the YAML file, and write it to the output directory.
+"""
+from ruamel.yaml import YAML
+import os
+import hashlib
+from pathlib import Path
+import sys
+
+# Read the AMMap configuration YAML file in a safe way
+yaml = YAML(typ='safe')
+with open(sys.argv[1], 'r') as f:
+    data = yaml.load(f)
+
+def constructCallable(name: str, hash: str, data: dict):
+    """Constructs the callable based on the name, hash, and data.
+
+    Args:
+        name: The name of the callable extracted from the YAML file.
+        hash: The truncated hash of the YAML file. Uses the default hashlib hashing algorithm (SHA256).
+        data: The data from the YAML file needed to set up the callable, including the elemental space components
+            and constraint min / max values to be enforced when evaluating the feasibility of a point.
+    """
+    # Extract the elements from the data and verify compliance
+    elements = data['elements']
+    assert isinstance(elements, list), 'Elements key must be a list'
+    assert all(isinstance(element, str) for element in elements), 'Elements key must be a list of strings'
+    assert len(elements) > 0, 'Elements key must have at least one element'
+
+    # Construct the output directory name
+    output = f"{name}_{hash}_{''.join(elements)}"
+    print(f"Constructing the callable: {output}")
+    
+    # Check if the directory exists, if not create it
+    if not os.path.exists(f'ammap/callables/{output}'):
+        os.makedirs(f'ammap/callables/{output}')
+        Path(f'ammap/callables/{output}/__init__.py').touch()
+
+    # Read the template file
+    with open('ammap/templates/EqScheil.py', 'r') as f:
+        template = f.read()
+
+    # Prepend the template with the constants based on the YAML file
+    constantsPayload = "# Constants based on the YAML file\n"
+    if 'elements' not in data:
+        raise ValueError('No elements key found in the YAML file')
+    else:
+        constantsPayload += f"ELEMENTS = {elements}\n"
+    
+    headerPayload = f'"""Linear Combination (atomic-fraction based) of elemental densities in {"-".join(elements)} system with constraints: '
+    
+    if 'min' in data:
+        constantsPayload += f"MIN = {data['min']}\n"
+        headerPayload += f"{data['min']}g/cm^3 < d"
+    if 'max' in data:
+        constantsPayload += f"MAX = {data['max']}\n"
+        if 'min' in data:
+            headerPayload += ' and '
+        headerPayload += f"d < {data['max']}g/cm^3"
+
+    headerPayload += '"""\n'
+    
+    payload = headerPayload + constantsPayload + '\n' + template
+
+    # Create the output file and write the payload
+    with open(f'ammap/callables/{output}/LCdensity.py', 'w') as f:
+        f.write(payload)
+
+# Check if the constraints key is present in the YAML file. Exit happily if not found.
+if 'constraints' not in data:
+    print('No constraints key found in the YAML file')
+    sys.exit(0)
+
+# Check if the LC density constraint is present in the YAML file. Exit happily if not found but raise an error if more than one is found.
+constraintMatch = [constraint['type'].lower().replace(' ', '') == 'lcdensity' for constraint in data['constraints']]
+if not any(constraintMatch):
+    print('No LC density constraint found in the YAML file')
+    sys.exit(0)
+elif sum(constraintMatch) > 1:
+    raise ValueError('More than one LC density constraint found in the YAML file')
+else:
+    pass
+
+# Check if the `elementalSpaces` list is present, is a list, and has at least one element. Exit with an error if not.
+if 'elementalSpaces' not in data:
+    raise ValueError('No elementalSpaces key found in the YAML file')
+elif not isinstance(data['elementalSpaces'], list):
+    raise ValueError('elementalSpaces key must be a list')
+elif len(data['elementalSpaces']) == 0:
+    raise ValueError('elementalSpaces key must have at least one element')
+
+# For each item under the `elementalSpaces` key, extract key information, elements, and construct the callable.
+constraintIndex = constraintMatch.index(True)
+constraint = data['constraints'][constraintIndex]
+name = data['name']
+hash = hashlib.sha256(str(data).encode()).hexdigest()[:6]
+
+for es in data['elementalSpaces']:
+    data = {'elements': es['elements']}
+    if 'min' in constraint:
+        data.update({'min': constraint['min']})
+    if 'max' in constraint:
+        data.update({'max': constraint['max']})
+    constructCallable(name, hash, data)
+
+
+        

ammap/callableBuilders/EqScheil.py

@@ -200,7 +200,7 @@ def getCD(temperature, solidFraction, CDPoints = [0.7,0.98], numDataThreshold =
 #        else:
 #            CD1.append(None)
 #            CD2.append(None)
-    return CD1, CD2
+#    return CD1, CD2
 
 def getNeighborCSC(temperature, solidFraction, CSCPoints=[0.4,0.9,0.99],numDataThreshold = 10):
     """Calculate the Critical Solidification Criteria.

ammap/callables/cracking.py

@@ -0,0 +1,105 @@
+from pycalphad import Database, equilibrium, variables as v
+from pycalphad.core.utils import instantiate_models, filter_phases, unpack_components
+from pycalphad.codegen.callables import build_phase_records
+import pandas as pd
+import math
+
+# Problem Specific setup for our 9-element space exploration. Make sure that the elements 
+# are in the same order as the composition vector that will be passed to the equilibrium_callable.
+dbf = Database("ammap/databases/TDB") #GLOBAL from yaml
+T = 'TEMPERATURE' #GLOBAL from yaml
+elementalSpaceComponents = 'ELEMENTS' #GLOBAL from yaml
+
+# Setup the pycalphad models
+phases = list(set(dbf.phases.keys()))
+comps = [s.upper() for s in elementalSpaceComponents]+['VA']
+phases_filtered = filter_phases(dbf, unpack_components(dbf, comps), phases)
+models = instantiate_models(dbf, comps, phases_filtered)
+phase_records = build_phase_records(dbf, comps, phases_filtered, {v.N, v.P, v.T}, models=models)
+expected_conds=[v.T]+[v.X(el) for el in comps[:-2]]
+default_conds={v.P: 101325, v.N: 1.0}
+
+# A neat callable for the equilibrium calculation that we will pass to the parallel graph exploration
+def equilibrium_callable(elP):
+    # Round to 6 decimal places, but make sure that 0.0 is not rounded to 0.0. pyCalphad will not like 
+    # if the components are exactly 0 or sum to exactly 1 in the equilibrium calculation for reasons
+    # that are beyond the scope of this tutorial.
+    elP_round = [round(v-0.000001, 6) if v>0.000001 else 0.0000001 for v in elP]
+    conds = {**default_conds, **dict(zip(expected_conds, [T] + elP_round))}
+    eq_res = equilibrium(
+        dbf, comps, phases_filtered, 
+        conds, model=models, phase_records=phase_records, calc_opts=dict(pdens=2000))
+    # Process the result into what we want -> a list of phases present
+    nPhases = eq_res.Phase.data.shape[-1]
+    phaseList = list(eq_res.Phase.data.reshape([nPhases]))
+    phasePresentList = [pn for pn in phaseList if pn!='']
+    if len(phasePresentList)==0:
+        # Decide on what to do if the equilibrium calculation failed to converge. By default, we will
+        # just pass and let the graph exploration scheme decide what to do about it.
+        # print(f"Point: {elP_round} failed to converge.")
+        pass
+    allphase = list(eq_res.Phase.data)
+    pFrac=eq_res.NP.data.shape[-1]
+    pFracList=list(eq_res.NP.data.reshape([nPhases]))
+    pFracPresent= [pn for pn in pFracList if not math.isnan(pn)]
+    return{
+    'Phases':phasePresentList,
+#    'allPhase': allphase,
+    'PhaseFraction':pFracPresent
+    }
+
+# Some extra code for the future tutorial on Scheil solidification :)
+from scheil import simulate_scheil_solidification
+# Meta Settings
+scheil_start_temperature = 'startTemperature' #GLOBAL from yaml
+liquid_phase_name = 'LIQUIDPHASE' #GLOBAL from yaml
+def scheil_callable(elP):
+    elP_round = [round(v-0.000001, 6) if v>0.000001 else 0.0000001 for v in elP]
+    initial_composition = {**dict(zip([v.X(el) for el in comps[:-2]], elP_round))}
+    
+    sol_res = simulate_scheil_solidification(
+        dbf, comps, phases_filtered, 
+        initial_composition, scheil_start_temperature, step_temperature=1.0)
+
+    phaseFractions = {}
+    for phase, ammounts in sol_res.cum_phase_amounts.items():
+        finalAmmount = round(ammounts[-1], 6)
+        if finalAmmount > 0:
+            phaseFractions[phase] = finalAmmount
+    test = sol_res.cum_phase_amounts
+    Lfrac = sol_res.fraction_liquid
+    Sfrac = sol_res.fraction_solid
+    scheilT = sol_res.temperatures
+    solT = sol_res.temperatures[len(scheilT) - 1]
+    ddict = sol_res.x_phases
+    keys_to_remove_from_ddict = []
+    for ddict_key, dict_value in ddict.items():
+        keys_to_remove_from_dict = [key for key, value in dict_value.items() if pd.isna(value).all()]
+    # Remove marked keys from the dictionary
+        for key in keys_to_remove_from_dict:
+            del dict_value[key]   
+    # If the dictionary is empty, mark the key in the defaultdict for removal
+        if not dict_value:
+            keys_to_remove_from_ddict.append(ddict_key)
+# Remove the marked keys from the defaultdict
+    for key in keys_to_remove_from_ddict:
+        del ddict[key]
+    yPhase = sol_res.Y_phases
+    for i, value in enumerate(Lfrac):
+        if value < 1:
+            break
+        liqT = sol_res.temperatures[i-1]
+    return {
+        'scheilT': scheilT,
+        'finalPhase': phaseFractions,
+        'Lfrac': Lfrac,
+        'Sfrac': Sfrac,
+        'solT': solT,
+        'liqT': liqT,
+        'phaseFractions': test,
+        'xPhase': ddict#,
+        #'yPhase': yPhase
+    }
+
+if __name__ == "main":
+    pass

ammap/templates/EqScheil.py

@@ -18,6 +18,16 @@ constraints:
     #  - SIGMA
   - type: LC density
     max: 11.5
+  - type: scheil
+    startTemperature: 2500
+    liquidPhase: 'LIQUID'
+  - type: cracking
+    criteria:
+     - FR
+     - Kou
+     - CSC
+     - iCSC
+     - sRDG
 
 elementalSpaces:
   - name: SS_V