ternary_frap.m

% Numerical solution of ternary FRAP model with solvent, bleached and
% unbleached species. Model is assumed to be equilibrated
% (bleached+unbleached=const.=pt). Then bleached species initial
% conditions are introduced. Integration of model via pdepe.

a = -50;
b = 0.025;
c = 1/10000;
u0 = 0.05;
e = 0.4;

%% Make useful mesh (by inverting the tanh profile and using this as spacing)
x = linspace(49.0, 51, 3000);
g = gamma0(x, a, 1*b, e);
g_unique = unique(g);
x = linspace(g_unique(1), g_unique(end-1), 300);
g_inv = spacing(x, a, 2*b, e);
g_inv = g_inv(2:end-1);
x = [linspace(49.0, g_inv(2), 30), g_inv(3:end-2), ...
    linspace(g_inv(end-1), 51, 30), linspace(51.1, 300, 300)];
%% Solve pde
tic
% x = [linspace(49.0, 49.8, 30), linspace(49.81, 50.2, 1000), ...
%      linspace(50.21, 51, 30) linspace(51.01, 200, 600)];
t = linspace(0.001, 1, 100);
fh_ic = @(x) flory_ic(x, a, u0);
fh_bc = @(xl, ul, xr, ur, t) flory_bc(xl, ul, xr, ur, t, u0);
fh_pde = @(x, t, u, dudx) flory_hugg_pde(x, t, u, dudx, a, b, e, c, u0);
sol = pdepe(2, fh_pde, fh_ic, fh_bc, x, t);
toc
%% Plotting
figure(1); hold on;
for i = 1:length(t)
    cla; xlim([49, 53]); ylim([0, 0.7]); 
    ax = gca;
    ax.FontSize = 16;
    xlabel('position'); ylabel('volume fraction');
    plot(x, phi_tot(x, a, b, e, u0), 'LineWidth', 2, 'LineStyle', '--'); 
    plot(x, sol(i, :), 'LineWidth', 2); pause(0.06);
% %     print([num2str(i),'.png'],'-dpng')
end

%% Plot and check derivatives of pt
figure; hold on;
x = linspace(40, 60, 100000);
plot(x, phi_tot(x, a, b, e, u0));
plot(x, gradient_analytical(x, a, b, e));
plot(x(1:end-1)+mean(diff(x))/2, ...
     diff(phi_tot(x, a, b, e, u0)/mean(diff(x))));
plot(x, gamma0(x, a, b, e));
figure; plot(gamma0(x, a, b, e), spacing(gamma0(x, a, b, e), a, b, e));
%% Function definitions for pde solver
function [c, f ,s] = flory_hugg_pde(x, t, u, dudx, a, b, e, c_p, u0)  
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot

pt = phi_tot(x, a, b, e, u0);
gra_a = gradient_analytical(x, a, b, e);
g0 = gamma0(x, a, b, e);
c = c_p;
f = g0*(1-pt)/pt*(pt*dudx-u*gra_a);
s = 0;
end

function u = flory_ic(x, a, u0)
    if x < -a
        u = 0.0;
    else
        u = u0;
    end
%     u0 = 0.3;
%  u0 = phi_tot(x, -50, 1);
end

function [pl,ql,pr,qr] = flory_bc(xl, ul, xr, ur, t, u0)
    pl = 0;
    ql = 1;
%     % No flux
%     pr = 0;%ur - 0.01;
%     qr = 1;
    % Dirichlet BC
    pr = ur - u0;
    qr = 0;
end

function g0 = gamma0(x, a, b, e)
    g0 = 1*e*(tanh((x+a)/b)+1)/2+0.001;
end

function sp = spacing(x, a, b, e)
    sp = b*atanh(2/(1*e)*(x-0.001)-1)-a;
end

function p = phi_tot(x, a, b, e, u0)
    p = e*(tanh(-(x+a)/b)+1)/2+u0;
end

function grad = gradient_analytical(x, a, b, e)
    grad = -e*(1-tanh(-(x+a)/b).^2)/b*0.5;
end