kawin.precipitation

Base class for precipitation models

The iteration method here may seem a bit odd, but it's mainly to reduce redundant calculations
For a model, we take that dX/dt = f(t, x), so only the time and current state is needed to compute dX/dt
While this is true for the KWN model, we need to perform quite a bit of calculations from x before being able to compute dX/dt
    This includes mass balance, nucleation rate and growth rate
For Euler:
    dX/dt is computed using the last state of the iteration, so in _calcDependentTerms, we just copy the last state in pData
For RK4:
    The first dX/dt is computed the same as for Euler
    Then for subsequent steps, we compute the new state using the current values of x (computed for each step of the iteration)
In both cases, after updating x, we compute the current state and append it to pData, which is then used for the next iteration

Parameters
----------
phases : list (optional)
    Precipitate phases (array of str)
    If only one phase is considered, the default is ['beta']
elements : list (optional)
    Solute elements in system
    Note: order of elements must correspond to order of elements set in Thermodynamics module
            Also, the list here should just be the solutes while the Thermodynamics module needs also the parent element
    If binary system, then default is ['solute']

Returns index of phase in list

Parameters
----------
phase : str (optional)
    Precipitate phase (defaults to None, which will return 0)

Resets simulation results
This does not reset the model parameters, however, it will clear any stopping conditions

Resets and initializes arrays for all variables
    time, temperature
    matrix composition, equilibrium composition (alpha and beta)
    driving force, impingement factor, nucleation barrier, critical radius, nucleation radius
    nucleation rate, precipitate density
    average radius, average aspect ratio, volume fraction

Extra variables include incubation offset and incubation sum

Time dependent variables will be set up as either
    (iterations)                     time, temperature
    (iterations, elements)           composition
    (iterations, phases, elements)   eq composition, total precipitate composition
    (iterations, phases)             Everything else
    This is intended for appending arrays to always be on the first axis

Converts precipitation data to dictionary

Converts dictionary of data to precipitation data

Appends new values to the variable list
NOTE: newVals must correspond to the same order as _packArrays with first axis as 1
    Ex rCrit is (n, phases) so corresponding new value should be (1, phases)
Since np append creates a new variable in memory, we have to reassign each term, then pack them into varList again
    TODO: it would be nice to reduce the number of times it copies, perhaps by preallocating some amount (say 1000)
            for each array and if we have not reached the end of the array, just stick the values at the latest index
            but once we reach the end of the array, we would append another 1000
            The after solving, we could clean up the arrays, or just use self.n to state where the end of the simulation is
I suppose we could make a list of str for each variable and call setattr

Sets constraints

Possible constraints:
---------------------
minRadius - minimum radius to be considered a precipitate (1e-10 m)
maxTempChange - maximum temperature change before lookup table is updated (only for Euler in binary case) (1 K)

maxDTFraction - maximum time increment allowed as a fraction of total simulation time (0.1)
minDTFraction - minimum time increment allowed as a fraction of total simulation time (1e-5)

checkTemperature - checks max temperature change (True)
maxNonIsothermalDT - maximum time step when temperature is changing (1 second)

checkPSD - checks maximum growth rate for particle size distribution (True)
maxDissolution - maximum relative volume fraction of precipitates allowed to dissolve in a single time step (0.01)

checkRcrit - checks maximum change in critical radius (False)
maxRcritChange - maximum change in critical radius (as a fraction) per single time step (0.01)

checkNucleation - checks maximum change in nucleation rate (True)
maxNucleationRateChange - maximum change in nucleation rate (on log scale) per single time step (0.5)
minNucleationRate - minimum nucleation rate to be considered for checking time intervals (1e-5)

checkVolumePre - estimates maximum volume change (True)
checkVolumePost - checks maximum calculated volume change (True)
maxVolumeChange - maximum absolute value that volume fraction can change per single time step (0.001)

minNucleateDensity - minimum nucleate density to consider nucleation to have occurred (1e-5)
dtScale - scaling factor to attempt to progressively increase dt over time

Sets function for beta calculation in binary systems
    1 for implementation seen in Perez et al, 2008 (default)
    2 for implementation similar to multicomponent systems

If using a multicomponent system, the beta function defaults to the 2nd
    So this function will not do anything

Parameters
----------
functionType : int
    ID for function
        1 for implementation seen in Perez et al, 2008 (default)
        2 for implementation similar to multicomponent systems

Returns Gibbs Thomson contribution of a particle given its radius

Parameters
----------
radius : float or array
    Precipitate radius
phase : str (optional)
    Phase to consider (defaults to first precipitate in list)

Adds condition to stop simulation when condition is met

Parameters
----------
condition: PrecipitateStoppingCondition
mode: str
    'or' or 'and
    Conditions with 'or' will stop the simulation when at least one condition is met
    Conditions with 'and' will stop the simulation when all conditions are met

Clears all stopping conditions

Sets up hidden parameters before solving
    Nucleation site density
    Grain boundary factors
    Strain energy

Overloads printHeader from GenericModel to do nothing
since status displays the necessary outputs

Prints various terms at latest step

Will print:
    Model time, simulation time, temperature, matrix composition
    For each phase
        Phase name, precipitate density, volume fraction, avg radius and driving force

Store array for non-derivative terms (which is everything except for the PBM models)

We use these terms for the first step of the iterators (for Euler, this is all the steps)
    For RK4, these terms will be recalculated in dXdt

Gets all dependent terms (everything but PBM variables) that are needed to find dXdt

Steps:
    1. Mass balance
    2. Driving force - must be done after mass balance to get the current matrix composition
    3. Growth rate - must be done after driving force since dG is needed in multicomponent systems
    4. Nucleation rate
    5. Nucleate radius - must be done after nucleation rate since derived classes can change nucleation rate

For the first iteration, self._currY will be None from the preProcess function, in this case, we want
    to just grab the latest values to avoid double calculations

Gets dXdt as a list for each phase
For the eulerian implementation, this is dn_i/dt for the bins in PBM for each phase

This is partially kind of dumb to have getdXdt and _getdXdt, however:
    getdXdt is to be compatible with GenericModel, which requires only x and t
    _getdXdt is to account for current PrecipitationData and growth rate as defined in this model

1) Updates internal arrays with new values of t and x
2) Updates particle size distribution
3) Updates coupled models
4) Check stopping conditions
5) Return new values and whether to stop the model

Euler implementation of KWN model

Parameters
----------
phases : list (optional)
    Precipitate phases (array of str)
    If only one phase is considered, the default is ['beta']
elements : list (optional)
    Solute elements in system
    Note: order of elements must correspond to order of elements set in Thermodynamics module
    If binary system, then defualt is ['solute']

Resets and initializes arrays for all variables

In addition to PrecipitateBase, the equilibrium aspect ratio area and population balance models are created here

Resets model results

Sets population balance model parameters for each phase

Parameters
----------
cMin : float
    Minimum bin size
cMax : float
    Maximum bin size
bins : int
    Initial number of bins
minBins : int
    Minimum number of bins - will not be used if adaptive = False
maxBins : int
    Maximum number of bins - will not be used if adaptive = False
adaptive : bool
    Sets adaptive bin sizes - bins may still change upon nucleation
phase : str
    Phase to consider (will set all phases if phase = None or 'all')

Sets recording parameters for PSD of specified phase

Parameters
----------
record : bool (optional)
    Whether to record PSD, defaults to True
phase : str (optional)
    Precipitate phase to record for
    Defaults to 'all', which will apply to all precipitate phases

Saves recorded PSD in npz format

Parameters
----------
filename : str
    File name to save to
    Note: the phase name will be added to the filename if all phases are being saved
compressed : bool (optional)
    Whether to save in compressed npz format
    Defualts to True
phase : str (optional)
    Phase to save PSD for
    Defaults to 'all', which will save a file for each phase

Loads particle size distribution for specified phase

Parameters
----------
data : array
    Array of data containing precipitate sizes
phase : str (optional)
    Phase to consider (defaults to first precipitate in list)

Returns PSD bounds of given phase

Parameters
----------
phase : str (optional)
    Phase to consider (defaults to first precipitate in list)

Returns Gibbs Thomson contribution of a particle given its radius

Parameters
----------
radius : array (optional)
    Precipitate radaii (defaults to None, which will use boundaries
        of the size classes of the precipitate PSD)
phase : str (optional)
    Phase to consider (defaults to first precipitate in list)

Returns population balance model of given phase

Parameters
----------
phase: str (optional)
    Phase to consider (defaults to first precipitate in list)

This creates a lookup table mapping the particle size classes to the interfacial composition

Sets up additional variables in addition to PrecipitateBase

Sets up additional outputs, population balance models, equilibrium aspect ratio and equilibrium compositions

Linear interpolation between self.eqAspectRatio and self.PBM[p].PSDbounds

Parameters
----------
R : float
    Equivalent spherical radius
p : int
    Phase index

The following checks are made
    1) change in number of particles moving between bins
        This is controlled by the implementation in PopulationBalanceModel,
        but essentially limits the number of particles moving between bins
    2) change in nucleation rate
        Time will be proportional to the 1/log(previous nuc rate / new nuc rate)
    3) change in temperature
        Limits how fast temperature can change
    4) change in critical radius
        Proportional to a percent change in critical radius
    5) estimated change in volume fraction
        Estimates the change in volume fraction from the nucleation rate and nucleation radius

Quick check to make sure particles below the thresholds are 0
    RdrivingForceIndex - only for binary, where energy from the Gibbs-Thompson effect is high enough
        that the free energy of the precipitate is above the free energy surface of the matrix phase
        and equilibrium cannot be calculated
    minRadius - minimum radius to be considered a precipitate

The _calcNucleationRate function in KWNBase calculates the nucleation rate as the
    probability that a site can form a nucleate that will continue to grow

To convert this probability to an actual nucleation rate, we multiply by the amount
    of available nucleation sites

The number of available sites is determined by:
    Available sites = All sites - used up sites + sites on parent precipitates
    The used up sites depends on the type of nucleation
        Bulk and grain corners - used sites = number of current precipitates
        Dislocation and grain edges - number of sites filled along the edges (assumes average radius of precipitates)
        Grain boundaries - number of sites filled along the faces (assumes average cross sectional area of precipitates)

Mass balance to find matrix composition with new particle size distribution
This also includes: volume fraction, precipitate density, average radius, average aspect ratio and sum of precipitate composition

Notes on computing composition from mass balance
    Concentration of the precipitates - needed to get matrix composition

    For a line compound with composition x^beta, this boils down to:
    x_0 = (1-f_v) * x^inf + f_v * x^beta
        Where x_0 is initial composition, f_v is volume fraction and x^inf is matrix composition

    For non-stoichiometric compounds, we want to integrate the precipitate composition as a function of radius
        We'll call this term f_conc (fraction + concentration of precipitates), so:
        x_0 = (1-f_v) * x^inf + f_conc
    
    For infinite precipitate diffusion, the concentration of a single precipitate is assumed to be homogenous
    f_conc = r_vol * vol_factor * sum(n_i * R_i^3 * x_i^beta)
        Where r_vol is V^alpha / V^beta and vol_factor is a factor for converting R^3 to volume (for sphere, this is 4*pi/3)

    For no diffusion in precipitate, the concentration depends on the history of the precipitate compositions and growth rate
    We just have to convert the summation to an integral of the time derivative of the terms inside
    f_conc = r_vol * vol_factor * sum(int(d(n_i * R_i^3 * x_i^beta)/dt, dt))
    We'll assume x_i^beta is constant with time (for 3 or more components, this is not true, but assume it doesn't change significantly per iteration - it'll also be a lot harder to account for)
    d(f_conc)/dt = r_vol * vol_factor * sum(d(n_i)/dt * R_i^3 * x_i^beta + 3 * R_i^2 * d(R_i)/dt * n_i * x_i^beta)
        d(n_i)/dt is the change in precipitates, since we don't record this, this is just (x[p] - self.PBM[p].PSD) / dt - with x[p] being the new number density for phase p
        d(R_i)/dt is the growth rate, however, since we use a eulerian solver, this corresponds to the growth rate of the bins themselves, which is 0
            If we were to use a langrangian solver, then d(n_i)/dt would be 0 (since the density in each bin would be constant) and d(R_i)/dt would be the growth rate at R_i
    Then we can calculate f_conc per iteration as a time integration like we do with some of the other variables

Returns current value of time and X
In this case, X is the particle size distribution for each phase

Returns dn_i/dt for each PBM of each phase

Corrects dXdt with the newly found dt, this adjusts the fluxes at the ends of the PBM so that we don't get negative bins

Calculates growth rate for a single phase
This is separated from _growthRateBinary since it's used in _calculatePSD

Matrix/precipitate composition are not calculated here since it's
already calculated in _createLookupBinary

Determines current growth rate of all particle size classes in a binary system

Calculates growth rate for a single phase
This is separated from _growthRateMulti since it's used in _calculatePSD

This will also calculate the matrix/precipitate composition 
for the radius in the PSD as well as equilibrium (infinite radius)

Determines current growth rate of all particle size classes in a multicomponent system

Updates particle size distribution with new x

Steps:
    1. Check if growth rate calculation failed with negative driving force
        We'll reset the PBM since we can't do much from here, but the chances of this happening should be pretty low
    2. Update the PBM with new x
    3. Check if the PBM needs to adjust the size class
        If so, then update the cached aspect ratio and precipitate composition with the new size classes
    4. Remove precipitates below a certain threshold (RdrivingForceIndex and minRadius)
    5. Calculate the dissolution index (index at which below are not considered when calculating dt)
        This is to prevent very small dt as the growth rate increases rapidly when R->0

Appends data from another PrecipitationData object to current one

Parameters for a single precipitate

Class for handling particle size distributions (PSD)
    This include time evolution, moments and graphing

NOTE: more formally, the PSD is defined such that the number of particles of size R in a volume
    is the integral of the PSD from R-1/2 to R+1/2
    For the discrete PSD, n_i = PSD_i * (R+1/2 - R-1/2)
    Here, we just store the PSD as the number of particles (n_i), but the formulation for
    growth rates and moments are the same since the factor (R+1/2 - R-1/2) cancels out

Parameters
----------
cMin : float (optional)
    Lower bound of PSD (defaults at 1e-10)
cMax : float (optional)
    Upper bound of PSD (defaults at 1e-9)
bins : int (optional)
    Initial number of bins (defaults at 150)
minBins : int (optional)
    Minimum number of bins over an order of magnitude (defaults at 100)
maxBins : int (optional)
    Maximum number of bins over an order of magnitude (defaults at 200)

Attributes
----------
originalMin : float
    Minimum bin size initially set
min : float
    Current minimum bin size (this will almost always be equal to originalMin)
originalMax : float
    Maximum bin size initially set
max : float
    Current maximum bin size
bins : int
    Default number of bins
minBins : int
    Minimum number of allowed bins
maxBins : int
    Maximum number of allowed bins
PSD : array
    Particle size distribution
PSDsize : array
    Average radius of each PSD size class
PSDbounds : array
    Radius at the bounds of each size class - length of array is len(PSD)+1

Resets the PSD to 0 and resets bin size and number of bins to original values
This will remove any size classes that were added since initialization

Enables recording of particle size distribution per iteration

The initial data in the recorded bin is t = 0, N_i = 0

The size of the recorded particle size distribution will be (n x max bins)
    Where n in the number of iterations
    max bins is the maximum number of bins, if the current number is smaller, the rest of the array will be 0

If recording, then reset the recorded bins to the original size (starting with t = 0, N_i = 0)
If not recording, then clear the recorded data

Disables recording

We won't clear the recorded bins here in case the user still wants to grab recorded data

Wrapper around enable and disable recording

Removes recorded data

Adds current PSD data to recorded arrays

TODO: Make sure this works when adaptive bins is False

Saves recorded data into npz format

Note: If recording is disabled, then this function will do nothing since
      there is nothing to save anyways

Parameters
----------
filename : str
    File name to save to

Loads recorded PSD

Returns PSD bounds, PSD bins and PSD from recorded data based off index

Since the number of bins is likely less than the max, we want to grab only the non-zero indices
TODO: two concerns
    1) this may remove the last 1 bins (this may be okay since we add new bins once the
        list bins has at least 1 particle), so the last bin would be 0 anyways

Sets particle size distribution to specific time if recorded

Parameter
---------
time : float
    Time to load PSD from, will load to nearest time available

For Euler implementation, sets whether to change the bin size when
the number of filled bins > maxBins or < minBins

If False, the bins will still be if nucleated particles are greater than the max bin size
and bins will still be added when the last bins starts to fill (but this will not change the bin size)

Sets constraints for minimum and maxinum number of bins over an order of magnitude

All bins will be overridden to an even number
Minimum number of bins will be overridden to be at most half of maximum bins
Default bins will be overridden to be halfway between min and max bins if out of range

Creates a particle size distribution from a set of data

Parameters
----------
data : array of floats
    Array of data to be inserted into PSD

Creates a particle size distribution from a function

Parameters
----------
function : function
    Takes in R and returns density

Stores current PSD and PSDbounds

Reverts to previous PSD and PSDbounds

NOTE: this appears to be unused
    (this was used in the previous KWNEuler implementation when the PSD could change within an iteration)
    (now it changes between iterations, so we don't need to revert back if something goes wrong)

Changes the size classes and resets the PSD

This is done by linear interpolation of the previous bins and PSD
And interpolating to the new bins and PSD
Due to differences in bin size (thus resolution of the PSD), the number density
    could be a little different. To correct for this, we get the 3rd moment of the
    previous PSD and the new PSD, and correct the new PSD to have the same 3rd moment

Parameters
----------
cMin : float
    Lower bound of PSD
cMax : float
    Upper bound of PSD
bins : int
    Number of bins
resetPSD : bool (optional)
    Whether to reset the PSD (defaults to False)

Adds an additional number of size classes to end of distribution

Parameters
----------
bins : int
    Number of bins to add

1) adds some bins to the end of the PSD if the last bin has at least 1 precipitate
    Number of bins is 1/4 of the original number of bins
2) If adaptive bin size is enabled, then two checks
    2a) if number of bins > max bins, then resize to have the number of bins be the minimum
    2b) if checking dissolution and number of filled bins < 1/2 min bins,
        then resize to last filled bin with the number of bins being the maximum

Parameters
----------
checkDissolution : bool
    Whether to check if the PSD is getting smaller and resize accordingly

Returns
-------
change : bool
    When the number of bins changed
newIndices : int or None
    The number of bins added to the PSD
    If the size of the bins changed, then this is None to indicate that resizing occured

Normalizes the PSD to 1

Normalizes the PSD with so that the moment of specified order will be 1

Parameters
----------
order : int (optional)
    Order of moment that PSD will be normalized to (defaults to 0)
    Using zeroth order will normalize the PSD so the sum of PSD will be 1

Finds indices when the volume fraction of particles below this index is
within the maximum amount (fraction-wise) that the PSD is allowed to dissolve

So find R_max where int(0, R_max, R^3 * dr) < maxDissolution * int(0, infinity, R^3 * dr)
The index is the correspoinding index to R_max

Parameters
----------
maxDissolution : float
    Max fraction allowed to dissolve
minIndex : int
    Minimum index which below, all particles are allowed to dissolve
    Upper limit on dissolution index

Returns
-------
max of [dissolution index, minIndex]

Calculates time interval for Euler implementation
    dt < dR / (2 * growth rate)
This ensures that at most, only half of particles in one size class can go to another

Also finds dt such that the max delta in growth rate is 0.4 dR
    We could use 0.5 dR which is the upper limit
        (for a given bin, the max change in density would remove all particles, with 0.5 getting smaller and 0.5 getting bigger)
    But 0.4 dR should be slightly more stable

TODO: allow variable ratio - this will make it more flexible for testing different time step constraints

Parameters
----------
currDT : float
    Current time interval, will be returned if it's smaller than what's given by the contraint
growth : array of floats
    Growth rate, must have lenth of bins+1 (or length of PSDbounds)
maxDissolution : float
    Maximum volume allowed to dissolve
startIndex : int
    First index to look at for growth rate, all indices below startIndex will be ignored
maxBinRatio : float (optional)
    Max ratio of particles in bin allowed to move to a nearby bin
    Default is 0.4

dn_i/dt = d(G*n)/dr + nucRate

d(G*n)/dr is calculated from two conditions
For positive growth rates - d(G*n)/dr|_i = n_i * flux_i / dr
For negative growth rates - d(G*n)/dr|_i = n_(i-1) * flux_i / dr
TODO : check that the two equations above represent the implementation

Parameters
----------
flux : numpy array (bins+1)
    Growth rate of particles in m/s
nucRate : float
    Nucleation rate in #/m^3/s
nucRadius : float
    Nucleation radius in m
psd : numpy array (bins)
    Particle size distribution with number density #/m3

Returns
-------
dXdt (bins) - corresponds to dn_i/dt

Given dt, correct the net flux so PSD will not be negative
    Essentially, the total number of particles leaving a bin should be less than or equal to the number of particles in the bin

Size of fluxes is bins+1 while for PSD, it is bins

For fluxes on the right (positive) side of the PSD
    We limit fluxes so that J_i+1 * dt < PSD_i

For fluxes on the left (negative) side of the PSD
    We limit fluxes so that -J_i * dt < PSD_i

Normally, this wouldn't be an issue since the time step from getDtEuler would limit J_i*dt to less than half the number of particles in a bin
    But because we set a dissolution threshold to ignore (to prevent extremely small dt since growth rate scales by 1/r), the bins below
    the dissolution threshold will not follow the constraint we apply in getDtEuler

Parameters
----------
dt : float
    time step
flux : numpy array (bins+1)
    Growth rate of particles in m/s
nucRate : float
    Nucleation rate in #/m^3/s
nucRadius : float
    Nucleation radius in m
psd : numpy array (bins)
    Particle size distribution with number density #/m3

Returns
-------
dXdt (bins) - corresponds to dn_i/dt corrected to avoid negative bins

Updates PBM with new values

Parameters
----------
time : float
    New time
newN : numpy array
    New number density

Given arbtrary PSD, return moment

Parameters
----------
order: float
    Moment order
weights: np.array (optional)
    Weights for each bin (must have length (bins,))
    If None, all bins are weighted equally
N: np.array (optional)
    PSD / number density (must have length (bins,))
    If None, then internal PSD will be used

Given arbtrary PSD, return cumulative moment (from 0 to max)

Parameters
----------
order: float
    Moment order
N: np.array (optional)
    PSD / number density (must have length (bins,))
    If None, then internal PSD will be used
weights: np.array (optional)
    Weights for each bin (must have length (bins,))
    If None, all bins are weighted equally

Sum of N

Length weighted moment of N

Area weighted moment of N

Volume weighted moment of N

Plots the N for each bin as a curve

Parameters
----------
pbm: PopulationBalanceModel
scale : float (optional)
    Scale factor for x-axis (defaults to 1)
    Note: this is for grain boundary nucleation where the
        reported precipitate radius differs from the radius
        determined by precipitate curvature
fill : bool (optional)
    Will fill area between PSD curve and x-axis (defaults to False)
ax: matplotlib Axes (optional)
*args, **kwargs - extra arguments for plotting

Plots the distribution density as a curve
Defined as f_i = N_i / (R_i+1/2 - R_i-1/2)
    This is such that N = \int(f dR)

Parameters
----------
pbm: PopulationBalanceModel
scale : float (optional)
    Scale factor for x-axis (defaults to 1)
    Note: this is for grain boundary nucleation where the
        reported precipitate radius differs from the radius
        determined by precipitate curvature
fill : bool (optional)
    Will fill area between PSD curve and x-axis (defaults to False)
ax: matplotlib Axes (optional)
*args, **kwargs - extra arguments for plotting

Plots cumulative size distribution density

Parameters
----------
axes : Axes
    Axis to plot on
scale : float (optional)
    Scale factor for x-axis (defaults to 1)
    Note: this is for grain boundary nucleation where the
        reported precipitate radius differs from the radius
        determined by precipitate curvature
order : int (optional)
    Moment of specified order
*args, **kwargs - extra arguments for plotting