kawin.precipitation.non_ideal

Defines class for calculating strain energy of a precipitate

Ellipsoidal precipitates will use the Eshelby's tensor
Spherical and Cuboidal precipitates will use the Khachaturyan's approximation

Sets method to invert the ohm term in calculating eshelby's tensor

Parameters
----------
method : str
    'numpy' - uses np.linalg.inv, which can be slower for batch, but runs through
                multiple checks for whether values are real/complex or if inverse exists
    'quick' - quick inverse using Cramer's rule assuming that values are real and inverse exists - recommended method

Sets resolution to which equilibrium aspect ratios are calculated

Equilibrium aspect ratios are found by calculated strain energy for a range of aspect ratios
and finding the aspect ratio giving the minimum energy (strain + interfacial energy)

If aspect ratio does not vary much in a given system, then the default parameters may lead to poor
prediction of the aspect ratios

Parameters
----------
resolution : float (optional)
    Minimum distance between aspect ratios when calculating cache (default at 0.01)
cachedRange : float (optional)
    Range of aspect ratio to calculate strain energy when updated cache (default at 5)

Sets method for calculating interfacial energy as a function of aspect ratio

Parameters
----------
method : str
    'eqradius' - interfacial energy is determined using the equivalent spherical radius
    'thermo' - interfacial energy is determined using dG/dSA (default)

Assumes spherical geometry for strain energy calculation
Uses Khachaturyan's approximation

Assumes cuboidal geometry for strain energy calculation
Uses Khachaturyan's approximation

Assumes ellipsoidal geometry for strain energy calculation
Uses Eshelby's tensor

If elastic strain energy is known to be a constant, this can be use to greatly
simplify calculations

Parameters
----------
energy - strain energy in J/m3

Sets elastic tensor of matrix using 2nd rank tensor

Parameters
----------
tensor : 6x6 array
    2nd rank elastic tensor

Sets elastic tensor of matrix by elastic constants, assuming cubic symmetry

Parameters
----------
c11 : float
    Modulus for compression
    c11 = E(1-nu) / (1+nu)(1-2nu)
c12 : float
    Modulus for dilation (accounts for compression and Poisson's ratio)
    c12 = E nu / (1+nu)(1-2nu)
c44 : float
    Modulus for shear
    c44 = (c11-c12)/2

Sets elastic tensor of matrix by 2 moduli

Parameters (only 2 has to be defined)
----------
E : float
    Elastic modulus
nu : float
    Poisson's ratio
G : float
    Shear modulus
lam : float
    Lame's first parameter
K : float
    Bulk modulus
M : float
    P-wave modulus

Sets elastic tensor of precipitate using 2nd rank tensor

Parameters
----------
tensor : 6x6 array
    2nd rank elastic tensor

Sets elastic tensor of precipitate by elastic constants, assuming cubic symmetry

Parameters
----------
c11 : float
    Modulus for compression
    c11 = E(1-nu) / (1+nu)(1-2nu)
c12 : float
    Modulus for dilation (accounts for compression and Poisson's ratio)
    c12 = E nu / (1+nu)(1-2nu)
c44 : float
    Modulus for shear
    c44 = (c11-c12)/2

Sets elastic tensor of precipitate by 2 moduli

Parameters (only 2 has to be defined)
----------
E : float
    Elastic modulus
nu : float
    Poisson's ratio
G : float
    Shear modulus
lam : float
    Lame's first parameter
K : float
    Bulk modulus
M : float
    P-wave modulus

Inputs 6x6 elastic matrix

Parameters
----------
tensor : matrix of floats
    Must be 6x6

Creates elastic tensor from c11, c12 and c44 constants assuming isotropic system

Parameters
----------
c11 : float
c12 : float
c44 : float

Set elastic tensor by defining at least two of the moduli
Everything will be converted to E, nu and G
If more than two parameters are defined, then the following priority will be used:
E - Youngs modulus
nu - Poisson's ratio
G - shear modulus
lam - Lame's first parameter
K - bulk modulus
M - P-wave modulus

NOTE: There's gotta be a better way to implement the conversions

Sets rotation matrix to be applied to the matrix

This is for cases where the axes of the precipitate does not align with the axes of the matrix
(e.g., the long/short axes of the precipitate is not parallel to the <100> directions of the matrix)

Parameters
----------
rot : matrix
    3x3 rotation matrix

Sets eigenstrain of precipitate

Parameters
----------
strain : float, array or matrix
    float - assume strain is the same along all 3 axis
    array - each index corresponds to strain in a given axis
    matrix - full 2nd rank strain tensor

NOTE: when using in conjunction with ShapeFactors, the axis are order specific
    For needle-like precipitates, x1, x2 and x3 correspond to (short axis, short axis, long axis)
    For plate-like precipitates, x1, x2 and x3 correspond to (long axis, long axis, short axis)

Sets applied stress tensor
Axes of stress tensor should be the same as the matrix

Parameters
----------
stress : float, array or matrix
    float - assume stress is the same along all 3 axis
    array - each index corresponds to stress in a given axis
    matrix - full 2nd rank stress tensor

NOTE: The applied stress is in reference to the coordinate axes of the precipitates
Thus, this is only valid if the following conditions are met:
    The matrix phase has a single orientation (either as a single crystal or highly textured)
    Precipitates will form only in a single orientation with respect to the matrix

TODO: It will be nice to expand on this.
    For polycrystalline matrix and randomly oriented precipitates, it should be possible
    to average the strain energy contributions over all matrix/precipitate orientations

Calculates applied strain tensor from applied stress tensor and elastic tensor

Sets up elastic constants

Converts 2nd rank elastic tensor to 4th rank

Parameters
----------
c : ndarray
    2nd rank elastic tensor

Returns
-------
c4 : ndarray
    4th rank elastic tensor

Converts 4th rank elastic tensor to 4nd rank

Parameters
----------
c4 : ndarray
    4th rank elastic tensor

Returns
-------
c : ndarray
    2nd rank elastic tensor

Inverts 4th rank tensor to give stiffness tensor

This is done by converting to 2nd rank, inverting, then converting back to 4th rank

Converts strain/stress vector to 2nd rank tensor

Parameters
----------
v : 1d array
    Strain/stress vector

Returns
-------
e : ndarray
    2nd rank elastic tensor

Converts 2nd rank tensor to vector

Parameter
---------
c : ndarray
    3x3 tensor

Returns
-------
v : 1darray
    Strain/stress vector

Rotates a 2nd rank tensor
T_ij = r_il * r_jk * T_lk

Parameters
----------
tensor : ndarray
    2nd rank tensor to rotate (3x3 array)

Rotates a 4th rank tensor
T_ijkl = r_im * r_jn * r_ok * r_lp * T_mnop

Parameters
----------
tensor : ndarray
    4th rank tensor to rotate (3x3x3x3 array)

Number of intervals to split domain along phi and theta for integration

Parameters
----------
phiInt : int
    Number of intervals to divide along phi
thetaInt : int
    Number of intervals to divide along theta
assumeSymmetric : bool (optional)
    If True (default), will only integrate Eshelby's tensor on a single quadrant and
    multiply the results by 8

Creates Lebedev quadrature points and nodes for integrating Eshebly's tensor
This is preferred over discretizing phi and theta

Parameters
----------
order : str
    'low' - uses quadrature order or 53 (974 points)
    'mid' - uses quadrature order or 83 (2354 points)
    'high' (default) - uses quadrature order or 131 (5810 points)

Unit normal vector of sphere

Parameters
----------
phi - azimuthal angle
theta - polar angle

Distance from center to surface of ellipsoid

Parameters
----------
a, b, c - radii of ellipsoid along x,y,z axes
phi, theta - azimuthal, polar angle

Ohm term for general system

Ohm_ij = inverse(C_iklj * n_k * n_l)

Hard coded inverse of m which is of shape (3,3,n)

numpy inv is more optimized for larger matrices, but can be slower for small
matrices such as a 2x2 or 3x3. We can take advantage of 3x3 matrices having a computable
inverse to make it faster

NOTE: this only works since we know that m has a shape of (3,3,n) and is only composed of real numbers

This function can probably be a bit more efficient, but quick 
profiling on sphInt gives around a 35x speedup compared to doing 
np.transpose(np.linalg.inv(np.transpose(m, (2,0,1))), (1,2,0)) 
where the slowdown was in np.linalg.inv

For matrix:
    |  a  b  c  |
    |  d  e  f  |
    |  g  h  i  |

Inverse is defined as:
    |  ei-fh  fg-di  dh-eg  |          |  A  B  C  |
    |  ch-bi  ai-cg  bg-ah  | / det -> |  D  E  F  | / det
    |  bf-ce  cd-af  ae-bd  |          |  G  H  I  |
    Where det = aA + bB + cC

Inverts ohm term using np.linalg.inv

numpy inverse function takes in an array of shape (m,n,n) and inverts each nxn matrix
    So we have to transpose m from (3,3,n) -> (n,3,3), then invert, then transpose (n,3,3) ->(3,3,n)

Faster version of calculating the Dijkl, which avoids
using a for loop to iterate across phi and theta

Dijkl term for Eshelby's theory

S_ijmn = -0.5 * C_lkmn * (D_iklj + D_jkli)

Multiplies 2 tensors
4th x 2nd -> c_ij = a_ijkl * b_kl
4th x 4th -> c_ijkl = a_ijmn * b_mnkl

u = -0.5 * V * sigma_ij * strain_ij

Strain energy of ellipsoidal particle with applied stress

Strain energy of ellipsoidal particle

Alternative method of strain energy on ellipsoidal particle using 2nd rank tensors

Strain energy of particle for when matrix and precipitate phases have different elastic tensors

Strain energy of particle for when matrix and precipitate phases have different elastic tensors using 2nd rank tensors

Khachaturyan's approximation for strain energy of spherical and cuboidal precipitates

Strain energy of perfect cube (cubic factor = sqrt(2))

Strain energy of perfect sphere (cubic factor = 1)

For cuboidal preciptitates, strain energy can be approximated
as a linear interpolation between a perfect cube and sphere

Generic strain energy function

Ohm term for cubic crystal symmetry

Needed for the Ohm term with cubic crystal symmetry

Needed for the Ohm term with cubic crystal symmetry

Equilibrium aspect ratio using golden ratio search

Parameters
----------
Rsph : float or array
    Equivalent spherical radius
gamma : float
    Interfacial energy
shpFactor : ShapeFactor object
a, b : float
    Min and max bounds
    Default is 1.001 and 100

Golden ratio search for a single radius

Update cached calculations

Clear cached calculations

Equilibrium aspect ratio by cached search

Parameters
----------
Rsph : float or array
    Equivalent spherical radius
gamma : float
    Interfacial energy
shpFactor : ShapeFactor object

Cached search for a single radius

Load lebedov nodes of specific order

Parameters
----------
order : int
    Can be 53, 83 or 131
    Determines the order of polynomial at which an exact integral solution can be found

Returns
-------
phi coordinates, theta coordinates, weights

Defines functions for shape factors of needle, plate and cuboidal precipitates
Shape factors involve the following:
    Aspect Ratio - ratio of long axis to short axis
    Eccentricity - characteristic value of the precipitate shape that depends on the aspect ratio
    Equivalent radius factor - correction factor to the radius to give the radius of a sphere with equivalent volume
    Kinetic factor - correction factor for the growth rate of the particle
    Thermodynamic factor / Area factor - correction factor for the critical driving force/radius for nucleation
    
    For a sphere (or needle and plate precipitates with aspect ratio of 1):
        Eccentricity = 0
        Equivalent radius factor = 1
        Kinetic factor = 1
        Thermodynamic factor = 1

NOTE:
normalRadaii, eqRadiusFactor, kineticFactor and thermoFactor 
    take in R and output the corresponding factors

_normalRadaiiEquation, __eqRadiusEquation, _kineticEquation and _thermoEquation
    take in aspect ratio and output the corresponding factors

Sets aspect ratio as a function of R (equivalent spherical radius)

Parameters
----------
ar : float or function
    Aspect ratio
    If function, must take in radius and return aspect ratio

Aspect ratio as a function of radius if the aspect ratio is set as a scalar

General shape setting function

Defualts to spherical

Sets factors for spherical precipitates

Factors for needle shaped precipitates

Factors for plate shaped precipitates

Factors for cuboidal shaped precipitates

Radius along the 3 axis to give a spherical volume of 1

Parameters
----------
R : float or array
    Equivalent spherical radius

Returns
-------
3 length array for each axis if R is scalar
n x 3 array if R is an array

Equivalent spherical radius factor for generic precipitates

Parameters
----------
R : float or array
    Equivalent spherical radius

Returns
-------
Eq. radius factor with same shape as R

Kinetic factor for generic precipitates

Parameters
----------
R : float or array
    Equivalent spherical radius

Returns
-------
Kinetic factor with same shape as R

Thermodynamic factor for generic precipitates

Parameters
----------
R : float or array
    Equivalent spherical radius

Returns
-------
Thermodynamic factor with same shape as R

Eccentricity given the aspect ratio (for needle and plate shaped precipitates)

Equivalent radius for a sphere (returns 1)

Equivalent radius for needle shaped precipitate

Equivalent radius for plate shaped precipitate

Equivalent radius for cuboidal shaped precipitate

Returns radius along the 3-axis for a volume of 1

Returns radius along the 3-axis for a volume of 1

Returns radius along the 3-axis for a volume of 1

Returns radius along the 3-axis for a volume of 1

Kinetic factor for a sphere (returns 1)

Kinetic factor for needle shaped precipitate

Kinetic factor for plate shaped precipitate

Kinetic factor for cuboidal shaped precipitate

Thermodynamic factor for a sphere (returns 1)

Thermodynamic factor for needle shaped precipitate

Thermodynamic factor for plate shaped precipitate

Thermodynamic factor for cuboidal shaped precipitate

Critical radius given a scalar aspect ratio

Critical radius given aspect ratio as a function of radius
Found by bisection method

Defines nucleation factors for bulk, dislocation, 
surface, edge and corner nucleation

This includes: surface area, volume, GB area removal, Rcrit and Gcrit

Attributes
----------
gbRemoval - factor to multiply by r**2 to get area of grain boundary removed
areaFactor - factor to multiply by r**2 to get surface area of precipitate
volumeFactor - factor to multiply by r**3 to get volume of precipitate

Sets nucleation site type

Parameters
----------
site - int
    Index for site type based off list on top of this file

Calculated area, volume and GB removal factors

Parameters
----------
gbEnergy - float
    Grain boundary energy
gamma - float
    Interfacial energy between precipitate and bulk

Critical radius for nucleation
This is only done for GB nucleation
    Bulk and dislocation nucleation needs to account for shape factors

R* = 2 * (b * \gamma_{\alpha \alpha} - a * \gamma_{\alpha \beta}) / (3 * c * dG)

Parameters
----------
dG - float
    Volumetric driving force

Critical driving force for nucleation

G* = 4/27 * (b * \gamma_{\alpha \alpha} - a * \gamma_{\alpha \beta})^3 / (c * dG)^2

or in terms of R*
G* = -R*^2 * (b * \gamma_{\alpha \alpha} - a * \gamma_{\alpha \beta} + c * r* * dG)

This is calculated in terms of R* since there is a check in the nucleation calculation that will
increase R* to minimum radius if R* is too small (very low nucleation barriers)

Parameters
----------
dG - float
    Volumetric driving force
Rcrit - float
    Critical radius for nucleation

Stores variables and functions to calculate effective diffusion distance

Sets up interpolation points for effective diffusion distance

Parameters
----------
n : int
    Number of points to interpolate between, two endpoints will be added
    at supersaturation of 0 and 1
lmin : float
    Minimum lambda value in log scale
lmax : float
    Maximum lambda value in log scale

Effective diffusion distance given supersaturation
This gives a constant to be multiplied by the particle radius

The effective diffusion distance is given by
eps = Q/2*lambda

Where Q is the super saturation and lambda is given by
Q = 2*lambda^2 - 2*lambda^3 * sqrt(pi) * exp(lambda^2) * (1 - erf(lambda))

Parameters
----------
supersaturation : float or array of floats
    (x - x_P) / (x_M - x_P) where x is matrix composition, x_M is parent composition and x_P is precipitate composition

Returns
-------
Effective diffusion distance (eps) to be used as (eps*R) in the growth rate equation

Lambda when Q approaches 0
This is done to prevent precision errors when multiplying exp*(1-erf)

Lambda when Q approaches 1
This is done to prevent precision errors when multiplying exp*(1-erf)

Effective diffusion distance given supersaturation
This gives a constant to be multiplied by the particle radius

When Q approaches 1, this equation oscillates due to floating point errors
Therefore, interpolation is performed between the two limits (Q -> 0 and Q -> 1)

This may be removed in later versions

Parameters
----------
supersaturation : float or array of floats
    (x - x_P) / (x_M - x_P) where x is matrix composition, x_M is parent composition and x_P is precipitate composition

Returns
-------
Effective diffusion distance (eps) to be used as (eps*R) in the growth rate equation

If neglecting the effective diffusion distance, return 1

Parameters
----------
supersaturation : float or array of floats
    (x - x_P) / (x_M - x_P) where x is matrix composition, x_M is parent composition and x_P is precipitate composition

Returns
-------
1 or array of 1s