%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File: pf3d_multi_aniso.m
%
%Purpose: 3-dimensional Multi Phase Field model (Steinbach 1996)
%         solving Ginzburg-Landau (or Cahn-Hillard) equation
%         by finite difference method
%         Anisotropic surface energy (cubic symmetry)
%        
% Wei Cai (caiwei@stanford.edu)
% Seunghwa Ryu (shryu@stanford.edu)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 3D multi phase field model with isotropic surface energy

if(~exist('N_grid'))
    N_grid = [50 50 50];
end

if(~exist('model_param'))
    %                h   M23 M31 M12 U23 U31 U12 E23 E31 E12 u1 u2 u3
    %                    Msv Mvl Mls Usv Uvl Uls Esv Evl Els ul us uv surface energy anisotropy
    model_param = [3 0.4 1.0 1.0 1.0 2.0 2.0 2.0 2.0 2.0 1.8 0  0  0  1.0 0.6993 -3.6595 -1.7384]; 
end

N_phases = model_param(1);    % 1: liquid, 2: solid, 3: vapor
h = model_param(2);           % grid size
M = zeros(N_phases,N_phases); % M_ik is mobility of interface i-k
U = zeros(N_phases,N_phases); % energy parameter of interface i-k
vareps = zeros(N_phases,N_phases); % energy parameter of interface i-k

M(2,3) = model_param(3); M(3,2) = model_param(3);
M(3,1) = model_param(4); M(1,3) = model_param(4);
M(1,2) = model_param(5); M(2,1) = model_param(5);

U(2,3) = model_param(6); U(3,2) = model_param(6);
U(3,1) = model_param(7); U(1,3) = model_param(7);
U(1,2) = model_param(8); U(2,1) = model_param(8);

vareps(2,3) = model_param(9);  vareps(3,2) = model_param(9);
vareps(3,1) = model_param(10); vareps(1,3) = model_param(10);
vareps(1,2) = model_param(11); vareps(2,1) = model_param(11);

mu = model_param(12:14);      % chemical potential of phases
m = zeros(N_phases,N_phases); % chemical potential difference between phases

eps0 = model_param(15);   % these modify Esv (implemented) and/or Els (not implemented)
eps1 = model_param(16);   % energy parameters
eps2 = model_param(17);   % energy parameters
eps3 = model_param(18);   % energy parameters

for i=1:3, for k=1:3, m(i,k) = mu(i) - mu(k); end; end

if(~exist('Niter'))
    Niter = 1000; % total number of simulation steps
end

if(~exist('dt'))
    dt = 5e-3; % time step
end

if(~exist('equation_type'))
    equation_type = 0; %0: Ginzburg-Landau, 1: Cahn-Hillard
end

if(~exist('plotfreq'))
    plotfreq = 50;
end

if(~exist('printfreq'))
    printfreq = 50;
end

if(~exist('solid_shape'))
    solid_shape = [1 1 1]; %sphere
    %solid_shape = [0 1 1]; %cylinder
end

if(~exist('solid_radius'))
    solid_radius = 0.25; 
end

cell_size = N_grid * h;
x = [0:N_grid(1)-1] * h;
y = [0:N_grid(2)-1] * h;
z = [0:N_grid(3)-1] * h;
Fdata = zeros(Niter,1);
Gdata = zeros(Niter,1);

if(~exist('phi'))
    phi{1} = zeros(N_grid); phi{2} = zeros(N_grid); phi{3} = zeros(N_grid);
    for i1=1:N_grid(1), for i2=1:N_grid, for i3=1:N_grid,
        if  (i1-N_grid(1)/2)^2*solid_shape(1) ...
           +(i2-N_grid(2)/2)^2*solid_shape(2) ...
           +(i3-N_grid(3)/2)^2*solid_shape(3) < (N_grid(1)*solid_radius)^2
            phi{2}(i1,i2,i3) = 1.0;
        end
    end, end, end
    for i1=1:N_grid(1), for i2=1:N_grid, for i3=1:N_grid,
       if  (i3<N_grid(3)*0.5)
           phi{1}(i1,i2,i3) = 1.0 - phi{2}(i1,i2,i3);
       end
    end, end, end
    phi{3} = 1 - phi{1} - phi{2};
end

pxind=[2:N_grid(1),1]; nxind=[N_grid(1),1:N_grid(1)-1];
pyind=[2:N_grid(2),1]; nyind=[N_grid(2),1:N_grid(2)-1];
pzind=[2:N_grid(3),1]; nzind=[N_grid(3),1:N_grid(3)-1];

%iteration begins
for iter = 1:Niter,
    
    for i=1:N_phases,
        dphidx{i} = (phi{i}(pxind,:,:)-phi{i}(nxind,:,:))/(2*h);
        dphidy{i} = (phi{i}(:,pyind,:)-phi{i}(:,nyind,:))/(2*h);
        dphidz{i} = (phi{i}(:,:,pzind)-phi{i}(:,:,nzind))/(2*h);
        d2phi{i}  = (phi{i}(pxind,:,:)+phi{i}(nxind,:,:) ...
                    +phi{i}(:,pyind,:)+phi{i}(:,nyind,:) ...
                    +phi{i}(:,:,pzind)+phi{i}(:,:,nzind)-6*phi{i})/(h^2);
        dphi_SQR{i} = (dphidx{i}.*dphidx{i}) + (dphidy{i}.*dphidy{i}) + (dphidz{i}.*dphidz{i});
        
        if i==2 % for solid phase only, compute interface orientation
            absdphi = sqrt(dphi_SQR{i} + 1e-24);
            dphi_CUBE = absdphi.^3;
            nx = dphidx{i} ./ absdphi;
            ny = dphidy{i} ./ absdphi;
            nz = dphidz{i} ./ absdphi;

            dnx_dphidx = 1./absdphi - dphidx{i} .* dphidx{i} ./ dphi_CUBE;
            dny_dphidx =            - dphidy{i} .* dphidx{i} ./ dphi_CUBE;
            dnz_dphidx =            - dphidz{i} .* dphidx{i} ./ dphi_CUBE;

            dnx_dphidy =            - dphidx{i} .* dphidy{i} ./ dphi_CUBE;
            dny_dphidy = 1./absdphi - dphidy{i} .* dphidy{i} ./ dphi_CUBE;
            dnz_dphidy =            - dphidz{i} .* dphidy{i} ./ dphi_CUBE;

            dnx_dphidz =            - dphidx{i} .* dphidz{i} ./ dphi_CUBE;
            dny_dphidz =            - dphidy{i} .* dphidz{i} ./ dphi_CUBE;
            dnz_dphidz = 1./absdphi - dphidz{i} .* dphidz{i} ./ dphi_CUBE;            

            triplesum = (nx.^2).*(ny.^2) + (ny.^2).*(nz.^2) + (nz.^2).*(nx.^2);
            EPS_loc   = eps0 + eps1.*triplesum + eps2.*(nx.*ny.*nz).^2 ...
                      + eps3*(triplesum.^2);
                  
            % the following is needed to compute gradients
            dEPS_dnx = 2.*eps1*nx.*((ny.^2)+(nz.^2)) ...
                     + 2.*eps2*nx.*(ny.^2).*(nz.^2)  ...
                     + 4.*eps3*nx.*((ny.^2)+(nz.^2)).*triplesum;
            dEPS_dny = 2.*eps1*ny.*((nz.^2)+(nx.^2)) ...
                     + 2.*eps2*ny.*(nz.^2).*(nx.^2)  ...
                     + 4.*eps3*ny.*((nz.^2)+(nx.^2)).*triplesum;
            dEPS_dnz = 2.*eps1*nz.*((nx.^2)+(ny.^2)) ...
                     + 2.*eps2*nz.*(nx.^2).*(ny.^2)  ...
                     + 4.*eps3*nz.*((nx.^2)+(ny.^2)).*triplesum;

            dEPS_dphidx = dEPS_dnx.*dnx_dphidx + dEPS_dny.*dny_dphidx + dEPS_dnz.*dnz_dphidx;
            dEPS_dphidy = dEPS_dnx.*dnx_dphidy + dEPS_dny.*dny_dphidy + dEPS_dnz.*dnz_dphidy;
            dEPS_dphidz = dEPS_dnx.*dnx_dphidz + dEPS_dny.*dny_dphidz + dEPS_dnz.*dnz_dphidz;
        end
    end
    
    F_pair = zeros(N_phases,N_phases); F = 0;
    for i=1:N_phases,
        for k=i+1:N_phases,
            grad_term = (phi{i}.*dphidx{k} - phi{k}.*dphidx{i}).^2 ...
                      + (phi{i}.*dphidy{k} - phi{k}.*dphidy{i}).^2 ...
                      + (phi{i}.*dphidz{k} - phi{k}.*dphidz{i}).^2 ;
            if i==2 && k==3 % only for solid-vapor interface
                EPS_ik = EPS_loc; % multiplication factor for anisotropy
            else 
                EPS_ik = 1.0;     % isotropic interface energy
            end
            F_pair(i,k) = vareps(i,k) * sum(sum(sum( 0.5 * EPS_ik.^2 .* grad_term ))) ...
                        + U(i,k) * sum(sum(sum( phi{i}.^2.*phi{k}.^2 )));
            F = F + F_pair(i,k);
        end
    end
    Fdata(iter)=F;
   
    % unmodified (original) rate of change for phi{i}
    for i=1:N_phases, dphidt0{i} = zeros(N_grid); end;
    for i=1:N_phases,
        for k=1:N_phases,
            if k ~= i
                if i==2 && k==3     % only for solid-vapor interface
                    dEPS_dphidx_loc = dEPS_dphidx;
                    dEPS_dphidy_loc = dEPS_dphidy;
                    dEPS_dphidz_loc = dEPS_dphidz;
                elseif i==3 && k==2 % only for solid-vapor interface (check with Seunghwa)
                    dEPS_dphidx_loc = -dEPS_dphidx;
                    dEPS_dphidy_loc = -dEPS_dphidy;
                    dEPS_dphidz_loc = -dEPS_dphidz;
                else
                    dEPS_dphidx_loc = 0;
                    dEPS_dphidy_loc = 0;
                    dEPS_dphidz_loc = 0;
                end
                grad_ik_x = phi{k}.*dphidx{i}-phi{i}.*dphidx{k};
                grad_ik_y = phi{k}.*dphidy{i}-phi{i}.*dphidy{k};
                grad_ik_z = phi{k}.*dphidz{i}-phi{i}.*dphidz{k};
                grad_ik_SQR = grad_ik_x.^2 + grad_ik_y.^2 + grad_ik_z.^2;
                tmp_x = (EPS_loc.^2).*grad_ik_x + grad_ik_SQR.*EPS_loc.*dEPS_dphidx_loc;
                tmp_y = (EPS_loc.^2).*grad_ik_y + grad_ik_SQR.*EPS_loc.*dEPS_dphidy_loc;
                tmp_z = (EPS_loc.^2).*grad_ik_z + grad_ik_SQR.*EPS_loc.*dEPS_dphidz_loc;    
                
                dFikdphii = 2*U(i,k)*phi{i}.*phi{k}.*(phi{k}-phi{i}) ...
                          + 4*m(i,k)*phi{i}.*phi{k} ...
                          - vareps(i,k) * (tmp_x(pxind,:,:)-tmp_x(nxind,:,:))/(2*h) ...
                          - vareps(i,k) * (tmp_y(:,pyind,:)-tmp_y(:,nyind,:))/(2*h) ...
                          - vareps(i,k) * (tmp_z(:,:,pzind)-tmp_z(:,:,nzind))/(2*h) ;
                dphidt0{i} = dphidt0{i} - M(i,k) * dFikdphii;
            end
        end
    end
    
    if (equation_type == 0)
        % Ginzburg-Landau
        for i=1:N_phases,
            dphidt{i} = dphidt0{i};
        end
    elseif (equation_type == 1)
        for i=1:N_phases,
            d2dphidt0{i} = (dphidt0{i}(pxind,:,:)+dphidt0{i}(nxind,:,:) ...
                           +dphidt0{i}(:,pyind,:)+dphidt0{i}(:,nyind,:) ...
                           +dphidt0{i}(:,:,pzind)+dphidt0{i}(:,:,nzind) ...
                           -6*dphidt0{i})/(h^2);
            dphidt{i} = d2dphidt0{i};
        end
    elseif (equation_type == 2)
        for i=1:N_phases,
            dphidt{i} = dphidt0{i} - mean(mean(mean(dphidt0{i})));
        end
    else
        disp(sprintf('Unknown equation_type %d',equation_type));
        return
    end
    
    for i=1:N_phases,
        phi{i} = phi{i} + dt* dphidt{i};
    end
    
    max_grad = zeros(N_phases,1);
    for i=1:N_phases,
        max_grad(i) = max(max(max(abs(dphidt{i}))));
    end
    Gdata(iter) = max(max_grad);
    
    if(mod(iter,printfreq)==0 || iter==1)
        disp(sprintf('iter=%d F=%11.3e G=%11.3e rel. vol. = (%.1f,%.1f,%.1f)%%',...
            iter,F,Gdata(iter),sum(sum(sum(phi{1})))/prod(N_grid)*100,...
            sum(sum(sum(phi{2})))/prod(N_grid)*100,...
            sum(sum(sum(phi{3})))/prod(N_grid)*100));
    end
    
    if(mod(iter,plotfreq)==0 || iter==1)
        figure(1);
        clf
        p_surf_1 = patch(isosurface(y,x,z,phi{1},0.5));
        isonormals(y,x,z,phi{1},p_surf_1)
        set(p_surf_1,'FaceColor','cyan','EdgeColor','none');
%         p_surf_2 = patch(isosurface(y,x,z,phi{2},0.5));
%         isonormals(y,x,z,phi{2},p_surf_2)
%         set(p_surf_2,'FaceColor','yellow','EdgeColor','none');
        p_surf_3 = patch(isosurface(y,x,z,phi{3},0.5));
        isonormals(y,x,z,phi{3},p_surf_3)
        set(p_surf_3,'FaceColor','magenta','EdgeColor','none');
        daspect([1 1 1])
        view(3); axis tight
        camlight 
        lighting gouraud
        xlim([min(y) max(y)]); ylim([min(x) max(x)]); zlim([min(z) max(z)]);
        alpha(0.7);
        
        figure(2);
        clf
        p_surf_2 = patch(isosurface(y,x,z,phi{2},0.5));
        isonormals(y,x,z,phi{2},p_surf_2)
        set(p_surf_2,'FaceColor','yellow','EdgeColor','none');
        daspect([1 1 1])
        view(3); axis tight
        camlight 
        lighting gouraud
        xlim([min(y) max(y)]); ylim([min(x) max(x)]); zlim([min(z) max(z)]);
        alpha(0.7);
        
        figure(3);
        subplot(2,1,1); 
        plot([0:iter-1]*dt,Fdata(1:iter));
        xlabel('t');
        ylabel('F');
        subplot(2,1,2);
        plot([0:iter-1]*dt,Gdata(1:iter));
        xlabel('t');
        ylabel('max(G)');        
    end
end