notebook.community

Edit and run 





In [1]:



    


TDB = """              
 ELEMENT /-   ELECTRON_GAS              0.0000E+00  0.0000E+00  0.0000E+00!
 ELEMENT VA   VACUUM                    0.0000E+00  0.0000E+00  0.0000E+00!
 ELEMENT CR   BCC_A2                    5.1996E+01  4.0500E+03  2.3560E+01!
 ELEMENT V    BCC_A2                    5.0941E+01  4.5070E+03  3.0890E+01!
 
 FUNCTION GHSERCR    2.98150E+02  -8856.94+157.48*T-26.908*T*LN(T)
     +.00189435*T**2-1.47721E-06*T**3+139250*T**(-1);  2.18000E+03  Y
      -34869.344+344.18*T-50*T*LN(T)-2.88526E+32*T**(-9);  6.00000E+03  N !

 FUNCTION GPCRBCC 298.14 0.0; 6000 N !

 FUNCTION GPCRLIQ 298.14 0.0; 6000 N !

 FUNCTION GCRFCC 298.14 +7284+.163*T+GHSERCR#; 6000 N !

 FUNCTION GHSERVV    2.98150E+02  -7930.43+133.346053*T-24.134*T*LN(T)
     -.003098*T**2+1.2175E-07*T**3+69460*T**(-1);  7.90000E+02  Y
      -7967.842+143.291093*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3;  
     2.18300E+03  Y
      -41689.864+321.140783*T-47.43*T*LN(T)+6.44389E+31*T**(-9);  
     4.00000E+03  N !

FUNCTION GLIQV      2.98150E+02  +12833.687+123.890501*T-24.134*T*LN(T)
     -.003098*T**2+1.2175E-07*T**3+69460*T**(-1)-5.19136E-22*T**7;  
     7.90000E+02  Y
      +12796.275+133.835541*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3
     -5.19136E-22*T**7;  2.18300E+03  Y
      -19617.51+311.055983*T-47.43*T*LN(T);  4.00000E+03  N !


FUNCTION UN_ASS 298.15 0; 300 N !
FUNCTION V001 298.15 0; 300 N !
FUNCTION V002 298.15 0; 300 N !
FUNCTION V003 298.15 0; 300 N !
FUNCTION V004 298.15 0; 300 N !
 

 TYPE_DEFINITION % SEQ *!
 DEFINE_SYSTEM_DEFAULT ELEMENT 2 !
 DEFAULT_COMMAND DEF_SYS_ELEMENT VA /- !


 PHASE LIQUID  %  1  1.0  !
   CONSTITUENT LIQUID  :CR,V :  !
   PARAMETER G(LIQUID,CR;0)  2.98150E+02  +15483.015+146.059775*T
    -26.908*T*LN(T)+.00189435*T**2-1.47721E-06*T**3+139250*T**(-1)
    +2.37615E-21*T**7;  2.18000E+03  Y
   -16459.984+335.616316*T-50*T*LN(T);  6.00000E+03  N REF1 !
   PARAMETER G(LIQUID,V;0)  2.98150E+02  +12833.687+123.890501*T
    -24.134*T*LN(T)-.003098*T**2+1.2175E-07*T**3+69460*T**(-1)
    -5.19136E-22*T**7;  7.90000E+02  Y
    +12796.275+133.835541*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3
    -5.19136E-22*T**7;  2.18300E+03  Y
    -19617.51+311.055983*T-47.43*T*LN(T);  4.00000E+03  N REF1 !
   PARAMETER L(LIQUID,CR,V;0)  2.98150E+02  -20224.86+V001#; 
    6.00000E+03  N REF3 !
   PARAMETER L(LIQUID,CR,V;1)  2.98150E+02  -13469.9+V002#; 
    6.00000E+03  N REF3 !


 PHASE BCC_A2  %&  2 1   3 !
   TYPE_DEFINITION & GES A_P_D BCC_A2 MAGNETIC  -1.0    4.00000E-01 !
   CONSTITUENT BCC_A2  :CR,V : VA :  !
   PARAMETER G(BCC_A2,CR:VA;0)  2.98150E+02  +GHSERCR#;   6.00000E+03   N 
    REF1 !
   PARAMETER TC(BCC_A2,CR:VA;0)  2.98150E+02  -311.5;
    6.00000E+03   N REF1 !
   PARAMETER BMAGN(BCC_A2,CR:VA;0)  2.98150E+02  -.008;   6.00000E+03   N 
    REF1 !
   PARAMETER L(BCC_A2,CR,V:VA;0) 298.14 -8253.85-3.61592*T+V003#;
    6000 N REF3 !
   PARAMETER L(BCC_A2,CR,V:VA;1) 298.14 +7494.82-8.69424*T+V004#;
    6000 N REF3 !  
   PARAMETER L(BCC_A2,CR,V:VA;2) 298.14 -17599.07+10.13142*T;
    6000 N REF3 !  
   PARAMETER L(BCC_A2,CR,V:VA;3) 298.14 +1042.58;
    6000 N REF3 !  
   PARAMETER G(BCC_A2,V:VA;0)  2.98150E+02  +GHSERVV#;
    4.00000E+03  N REF1 !

   
 PHASE HCP_A3  %)  2 1   .5 !
   TYPE_DEFINITION ) GES A_P_D HCP_A3 MAGNETIC  -3.0    2.80000E-01 !
   CONSTITUENT HCP_A3  :CR,V : VA :  !
   PARAMETER G(HCP_A3,CR:VA;0)  2.98150E+02  -4418.94+157.48*T
    -26.908*T*LN(T)+.00189435*T**2-1.47721E-06*T**3+139250*T**(-1);  
    2.18000E+03  Y
    -30431.344+344.18*T-50*T*LN(T)-2.88526E+32*T**(-9);
    6.00000E+03  N REF1 !
   PARAMETER TC(HCP_A3,CR:VA;0)  2.98150E+02  -1109; 
    6.00000E+03   N REF1 !
   PARAMETER BMAGN(HCP_A3,CR:VA;0)  2.98150E+02  -2.46;
    6.00000E+03   N REF1 !
   PARAMETER G(HCP_A3,V:VA;0)  2.98150E+02  -3930.43+135.746053*T
    -24.134*T*LN(T)-.003098*T**2+1.2175E-07*T**3+69460*T**(-1);
    7.90000E+02  Y
    -3967.842+145.691093*T-25.9*T*LN(T)+6.25E-05*T**2-6.8E-07*T**3;  
    2.18300E+03  Y
    -37689.864+323.540783*T-47.43*T*LN(T)+6.44389E+31*T**(-9);
    4.00000E+03  N REF1 !
   PARAMETER L(HCP_A3,CR,V:VA;0)  2.98150E+02  +5000; 
    6.00000E+03  N REF3 !


"""



    











In [2]:



    


from pycalphad import Database, Model, equilibrium
from pycalphad.core.solver import InteriorPointSolver
import pycalphad.variables as v

class ProblemSaver(InteriorPointSolver):
    saved_problem = [None]
    def solve(self, prob):
        self.saved_problem[0] = prob
        return super(ProblemSaver, self).solve(prob)



    











In [3]:



    


dbf = Database(TDB)
comps = ['CR', 'V', 'VA']
phases = list(dbf.phases.keys())
conds = {v.T: 2500, v.P: 101325, v.X('V'): 0.25}
equilibrium(dbf, comps, phases, conds, solver=ProblemSaver(),
            parameters={'V001': 0,'V002': 0,'V003': 0,'V004': 0})



    


















    
Out[3]:








<xarray.Dataset>
Dimensions:    (N: 1, P: 1, T: 1, X_V: 1, component: 2, internal_dof: 3, vertex: 3)
Coordinates:
  * N          (N) float64 1.0
  * P          (P) float64 1.013e+05
  * T          (T) float64 2.5e+03
  * X_V        (X_V) float64 0.25
  * vertex     (vertex) int64 0 1 2
  * component  (component) <U2 'CR' 'V'
Dimensions without coordinates: internal_dof
Data variables:
    NP         (N, P, T, X_V, vertex) float64 1.0 nan nan
    GM         (N, P, T, X_V) float64 -1.757e+05
    MU         (N, P, T, X_V, component) float64 -1.644e+05 -2.099e+05
    X          (N, P, T, X_V, vertex, component) float64 0.75 0.25 ... nan nan
    Y          (N, P, T, X_V, vertex, internal_dof) float64 0.75 0.25 ... nan
    Phase      (N, P, T, X_V, vertex) <U6 'LIQUID' '' ''
Attributes:
    engine:   pycalphad 0.8.1.post1+4.g198fd807
    created:  2020-01-31T22:07:12.356316

















In [4]:



    


import numpy as np
ProblemSaver.saved_problem[0].chemical_potential_parameter_gradient(ProblemSaver.saved_problem[0].x0)



    


















    
Out[4]:








array([[0.18674624, 0.09393448, 0.        , 0.        ],
       [0.18674624, 0.09393448, 0.        , 0.        ]])

















In [5]:



    


import numpy as np

def error_gradient(temperature, x_solid, x_liquid, parameters):
    param_values = np.array([value for key, value in sorted(parameters.items())])
    eq_solid_point = equilibrium(dbf, comps, ['BCC_A2', 'LIQUID'],
                                 {v.T: temperature, v.P: 101325, v.X('V'): x_solid},
                                 solver=ProblemSaver(), parameters=parameters)
    eq_solid_MU_grad = ProblemSaver.saved_problem[0].chemical_potential_parameter_gradient(ProblemSaver.saved_problem[0].x0)
    eq_liquid_point = equilibrium(dbf, comps, ['BCC_A2', 'LIQUID'],
                                 {v.T: temperature, v.P: 101325, v.X('V'): x_liquid},
                                 solver=ProblemSaver(), parameters=parameters)
    eq_liquid_MU_grad = ProblemSaver.saved_problem[0].chemical_potential_parameter_gradient(ProblemSaver.saved_problem[0].x0)
    candidate_hyperplane = np.nanmean([eq_solid_point.MU.values.squeeze(),
                                      eq_liquid_point.MU.values.squeeze()], axis=0)
    candidate_hyperplane_grad = np.nanmean([eq_solid_MU_grad,
                                            eq_liquid_MU_grad], axis=0)
    eq_single_phase_solid = equilibrium(dbf, comps, ['BCC_A2'],
                                 {v.T: temperature, v.P: 101325, v.X('V'): x_solid},
                                 solver=ProblemSaver(), parameters=parameters)
    eq_solid_single_phase_MU_grad = ProblemSaver.saved_problem[0].chemical_potential_parameter_gradient(ProblemSaver.saved_problem[0].x0)
    result_error = np.sum((eq_single_phase_solid.MU.values.squeeze() - candidate_hyperplane) * \
                        [1-x_solid, x_solid]) ** 2
    result_grad = 2 * np.sqrt(result_error) * np.dot((eq_solid_single_phase_MU_grad - candidate_hyperplane_grad).T,
                                                     [1-x_solid, x_solid])
    eq_single_phase_liquid = equilibrium(dbf, comps, ['LIQUID'],
                                 {v.T: temperature, v.P: 101325, v.X('V'): x_liquid},
                                 solver=ProblemSaver(), parameters=parameters)
    eq_liquid_single_phase_MU_grad = ProblemSaver.saved_problem[0].chemical_potential_parameter_gradient(ProblemSaver.saved_problem[0].x0)
    result_error_2 = np.sum((eq_single_phase_liquid.MU.values.squeeze() - candidate_hyperplane) * \
                         [1-x_liquid, x_liquid]) ** 2
    result_error += result_error_2
    result_grad += 2 * np.sqrt(result_error_2) * np.dot((eq_liquid_single_phase_MU_grad - candidate_hyperplane_grad).T,
                                                        [1-x_liquid, x_liquid])
    return result_error, result_grad



    











In [6]:



    


%matplotlib inline
import matplotlib.pyplot as plt

params_to_try = np.linspace(-10000, 10000, num=100)
result = list(error_gradient(2044, 0.3, 0.3, {'V001': x,'V002': 0,'V003': 0,'V004': 0}) for x in params_to_try)
errors = [x[0] for x in result]
error_grad = [x[1][0] for x in result]



    











In [7]:



    


plt.plot(params_to_try, errors)
plt.xlabel('V001')
plt.ylabel('Error')



    


















    
Out[7]:








Text(0, 0.5, 'Error')










    
























In [8]:



    


plt.plot(params_to_try, error_grad)
plt.xlabel('V001')
plt.ylabel('Error Gradient')



    


















    
Out[8]:








Text(0, 0.5, 'Error Gradient')










    
























In [ ]:

Content source: richardotis/pycalphad-sandbox

Similar notebooks:

notebook.community | gallery | about