verifyStabilization.m

function all_stable = verifyStabilization(sol_matrix, t_array, time_fraction)
%VERIFYSTABILIZATION - Checks if the solution reached a stabilized status comparing the final time point with a manually specified time point
% beware that this can have false positives when the analyzed solution had
% a tmax so small that the ionic movement did not even start
%
% Syntax:  all_stable = verifyStabilization(sol_matrix, t_array, time_fraction)
%
% Inputs:
%   SOL_MATRIX - a solution matrix, not the full struct, as contained in
%     structs created by PINDRIFT
%   T_ARRAY - the time mesh array, as generated by meshgen_t and included
%     in structs created by PINDRIFT
%   TIME_FRACTION - the fraction of the final time that can be used as a
%     reference. In case of a stabilization, the smaller the stricter the
%     check. In case of a periodic time evolution, a specific value
%     corresponding to the previous-to last cycle have to be used. The
%     value needs to be strictly between 0 and 1.
%
% Outputs:
%   ALL_STABLE - a logic asserting if the solution in input reached its
%     stability
%
% Example:
%   verifyStabilization(ssol_i_1S_SR.sol, ssol_i_1S_SR.t, 1e-3)
%     check if the solution did not change much from a time which is one
%     thousandth of the final time point and the final time point
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also pindrift.

% Author: Ilario Gelmetti, Ph.D. student, perovskite photovoltaics
% Institute of Chemical Research of Catalonia (ICIQ)
% Research Group Prof. Emilio Palomares
% email address: iochesonome@gmail.com
% Supervised by: Dr. Phil Calado, Dr. Piers Barnes, Prof. Jenny Nelson
% Imperial College London
% October 2017; Last revision: January 2018

%------------- BEGIN CODE --------------

% name of the variables
names = ["electrons", "holes", "ions", "potential"];
% which values have to be considered in a linear or in a log10 scale
compare_log = [true, true, false, false];
all_stable = true;

% no need to calculate end_time for each of the 4 solutions: if they
% break they break at the same time
% using t(end) could get a time out of the solution when the computation broke before reaching the final time
end_time = t_array(length(sol_matrix(:, 1, 1)));
[~, time_index] = min(abs(t_array - time_fraction * end_time)); % get time mesh index at a specified percentage of maximum time reached

% for each component of the solution verify that didn't change significantly over the requested time
for i = 1:length(sol_matrix(1, 1, :))
    profile_at_time = sol_matrix(time_index, :, i); % take profile of values at a certain time of evolution
    profile_end = sol_matrix(end, :, i); % take profile of values at the end of time
    if compare_log(i) % for variables ranging on huge scales comparing the log values makes more sense
        profile_at_time = log10(profile_at_time);
        profile_end = log10(profile_end);
    end

    difference = sum(abs(profile_end - profile_at_time)); % sum up all the differences between the profiles
    % if the values change more than one tenthousandth, the solution is
    % considered not stable
    threshold = 1e-4 * sum(abs(profile_end - mean(profile_end))); % sum up absolute values, ignore constant bias
    stable = difference <= threshold;

    if ~stable
        warning('pindrift:verifyStabilization',...
            'Comparing final solutions at %s s and %s s showed that the %s distribution did not reach stability. Consider trying with a greater tmax.',...
            num2str(t_array(time_index)), num2str(t_array(end)), names(i));
    end
    % true just if all the variables are stabilized
    all_stable = all_stable && stable;
end

if all_stable
    disp("Stabilisation verified");
end

%------------- END OF CODE --------------