Contents

Preallocating

Monte Carlo diffusion calculation

Preallocating

tic % track total time to track computational efficiency 
timeout = 5e5; % maximum number of time steps allowed for computational efficiency 
maxout = 0; % counter for parameter combinations that reached the maximum number of time steps without converging on a solution, begins at 0 

% arrays to hold best fitting results considering only fit for isotopic composition for each paramter combination 
Dbestd7 = zeros(it, length(x)); % isotope composition 
DbestC = zeros(it, length(x)); % concentration 
Dfintime = zeros(it, 1); % time 
Dbestfit = zeros(it, 1); % fit statistic 

% arrays to hold best fitting results considering concentration and isotopic composition for each paramter combination 
Tbestd7 = zeros(it, length(x)); % isotope composition 
TbestC = zeros(it, length(x)); % concentration 
Tfintime = zeros(it, 1); % time 
Tbestfit = zeros(it, 1); % fit statistic 

parfor_progress(it) % monitor of parallel for loop progress 

Monte Carlo diffusion calculation

parfor i = 1:it % beginning parallel for loop 
    parfor_progress; 
    % Temporary variables for the current iteration 
    C7 = zeros(timeout + 1, length(x)); % array for concentration of 7Li 
    C6 = zeros(timeout + 1, length(x)); % array for concentration of 6Li 

    % Initial concentration gradients 
    C7(1, :) = Ci7(i, 1) * (x < 0) + Ci7(i, 2) * (x > 0); % setting the initial rind (x < 0) and block (x > 0) concentrations 
    C6(1, :) = Ci6(i, 1) * (x < 0) + Ci6(i, 2) * (x > 0); % setting the initial rind (x < 0) and block (x > 0) concentrations 

    d7 = zeros(timeout, length(x)); % array for isotopic composition at each time step 
    C = zeros(timeout, length(x)); % array for total concentration at each time step 
    modfitd7 = zeros(timeout, 1); % array for fit statistic considering isotopic composition 
    modfitC = zeros(timeout, 1); % array for fit statistic considering concentration 

    % Numerical calculation for each time step 
    for n = 1:timeout 
        % Boundary conditions 
        C7(n + 1, x < 0) = C7(n, 2); % pinned boundary condition in the rind 
        % C7(n + 1, 1) = C7(n, 2); % open boundary condition used for the rind in the unpinned model 
        C7(n + 1, length(x)) = C7(n, length(x) - 1); % open boundary condition in the block 
        C6(n + 1, x < 0) = C6(n, 2); % pinned boundary condition in the rind 
        % C6(n + 1, 1) = C6(n, 2); % open boundary condition used for the rind in the unpinned model 
        C6(n + 1, length(x)) = C6(n, length(x) - 1); % open boundary condition in the block 

        % Numerical calculation for each spatial coordinate 
        C7(n + 1, 19:length(x) - 1) = C7(n, 19:length(x) - 1) + ((De7(i) * dt(i) ./ (Ke(i, 19:length(x) - 1) * dx^2)) .* (C7(n, 18:length(x) - 2) - 2 *
        C6(n + 1, 19:length(x) - 1) = C6(n, 19:length(x) - 1) + ((De6(i) * dt(i) ./ (Ke(i, 19:length(x) - 1) * dx^2)) .* (C6(n, 18:length(x) - 2) - 2 *

        % Convert isotope concentration gradients to d7Li and total C 
        d7(n, :) = ((((C7(n + 1, :) ./ C6(n + 1, :)) * (m6/m7)) ./ LSVEC76) - 1) * 1000; % ‰ relative to LSVEC for comparison with measured data 
        C(n, :) = (C7(n + 1, :) + C6(n + 1, :)) .* Kd; % total concentration = concentration of 7Li + concentration of 6Li 

        % Assess the model fit 
        % fit relative to measured isotopic data 
        modfitd7(n) = sum((((((d7(n, match) < measu) .* d7(n, match)) - ((d7(n, match) < measu) .* measu)).^2) ./ (meas(3, :).^2)) + (((((d7(n, match) 

        % fit relative to measured concentrations 
        modfitC(n) = sum((((((C(n, match) < measuC) .* C(n, match)) - ((C(n, match) < measuC) .* measuC)).^2) ./ (meas(2, :).^2)) + (((((C(n, match) > 

        % Check for good fit conditions 
        % all models run at least 2000 time steps 
        % if both isotopic and concentration fits are at least 5% worse than the best fit so far, the model is allowed to stop and the best fit is reta
        if n > 2000 && ((modfitd7(n) - min(modfitd7(1:n))) / min(modfitd7(1:n))) > 0.05 && ((modfitC(n) - min(modfitC(1:n))) / min(modfitC(1:n))) > 0.0
            % Store results 
            [~, I] = min(modfitd7(1:n), [], 1); 
            Dbestd7(i, :) = d7(I, :); 
            DbestC(i, :) = C(I, :); 
            Dfintime(i) = I * dt(i); 
            Dbestfit(i) = modfitd7(I); 

            [~, I] = min((modfitd7(1:n) + modfitC(1:n)), [], 1); 
            Tbestd7(i, :) = d7(I, :); 
            TbestC(i, :) = C(I, :); 
            Tfintime(i) = I * dt(i); 
            Tbestfit(i) = modfitd7(I) + modfitC(I); 
            break; 



       elseif C(n,size(C,2))>11.8 || C(n,1)<13.5 || d7(n,size(C,2))>-0.9 
            % If the concentration in the block or rind no longer falls within the measured range, or the isotopic composition in the block exceeds the
            % d7Li fit 
            [~, I] = min(modfitd7(1:n), [], 1); 
            Dbestd7(i, :) = d7(I, :); 
            DbestC(i, :) = C(I, :); 
            Dfintime(i) = I * dt(i); 
            Dbestfit(i) = modfitd7(I); 

            % total fit 
            [~, I] = min((modfitd7(1:n) + modfitC(1:n)), [], 1); 
            Tbestd7(i, :) = d7(I, :); 
            TbestC(i, :) = C(I, :); 
            Tfintime(i) = I * dt(i); 
            Tbestfit(i) = modfitd7(I) + modfitC(I); 
            break 
        elseif n==timeout 
            % if the number of time steps reaches the timeout value set of computational efficiency, the model is stopped and the best fit is retained 
            % d7Li fit 
            [~, I] = min(modfitd7(1:n), [], 1); 
            Dbestd7(i, :) = d7(I, :); 
            DbestC(i, :) = C(I, :); 
            Dfintime(i) = I * dt(i); 
            Dbestfit(i) = modfitd7(I); 

            % total fit 
            [~, I] = min((modfitd7(1:n) + modfitC(1:n)), [], 1); 
            Tbestd7(i, :) = d7(I, :); 
            TbestC(i, :) = C(I, :); 
            Tfintime(i) = I * dt(i); 
            Tbestfit(i) = modfitd7(I) + modfitC(I); 
            maxout=[maxout i]; 
        end 
    end 
end 
parfor_progress(0); 
save results_pin_A10.mat 
toc 

Published with MATLAB® R2024a