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

% 3D single phase field model with isotropic surface energy

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

if(~exist('model_param'))
    % surface energy parameters cooresponding to anisotropy of Si 
    model_param = [0.4 1.0 1.0 1.0 0.6993 -3.6595 -1.7384];
end

h = model_param(1);      % grid size
M = model_param(2);      % kinetic coefficient
U = model_param(3);      % energy parameters
eps0 = model_param(4);   % energy parameters
eps1 = model_param(5);   % energy parameters
eps2 = model_param(6);   % energy parameters
eps3 = model_param(7);   % energy parameters


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

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

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

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

if(~exist('printfreq'))
    printfreq = 50;
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 = zeros(N_grid);
    for i=1:N_grid(1), for j=1:N_grid(2), for k=1:N_grid(3),
       if (i-N_grid(1)/2)^2+(j-N_grid(2)/2)^2+(k-N_grid(3)/2)^2 < (N_grid(1)*0.3)^2
           phi(i,j,k) = 1.0;
       end
    end; end; end
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,
    dphidx = (phi(pxind,:,:)-phi(nxind,:,:))/(2*h);
    dphidy = (phi(:,pyind,:)-phi(:,nyind,:))/(2*h);
    dphidz = (phi(:,:,pzind)-phi(:,:,nzind))/(2*h);
    d2phi  = (phi(pxind,:,:)+phi(nxind,:,:)+phi(:,pyind,:)+phi(:,nyind,:) ...
             +phi(:,:,pzind)+phi(:,:,nzind)-6*phi)/(h^2);

    dphi_SQR = (dphidx.*dphidx) + (dphidy.*dphidy) + (dphidz.*dphidz);
    absdphi = sqrt(dphi_SQR + 1e-24);
    dphi_CUBE = absdphi.^3;
    
    nx = dphidx ./ absdphi;
    ny = dphidy ./ absdphi;
    nz = dphidz ./ absdphi;
    
    dnx_dphidx = 1./absdphi - dphidx .* dphidx ./ dphi_CUBE;
    dny_dphidx =            - dphidy .* dphidx ./ dphi_CUBE;
    dnz_dphidx =            - dphidz .* dphidx ./ dphi_CUBE;

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

    dnx_dphidz =            - dphidx .* dphidz ./ dphi_CUBE;
    dny_dphidz =            - dphidy .* dphidz ./ dphi_CUBE;
    dnz_dphidz = 1./absdphi - dphidz .* dphidz ./ 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);
                     
    % vareps in 'iso.m' now becomes (0.5 EPS_loc.^2)
    F = (U * sum(sum(sum( phi.^2.*(1-phi).^2 ))) ...
        + 0.5 * sum(sum(sum( EPS_loc.^2 .* dphi_SQR)))) * (h^3);
    Fdata(iter)=F;
   
    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;

    tmp_x = (EPS_loc.^2).*dphidx + dphi_SQR.*EPS_loc.*dEPS_dphidx;
    tmp_y = (EPS_loc.^2).*dphidy + dphi_SQR.*EPS_loc.*dEPS_dphidy;
    tmp_z = (EPS_loc.^2).*dphidz + dphi_SQR.*EPS_loc.*dEPS_dphidz;    
    
    dFdphi = 2*U*phi.*(1-phi).*(1-2*phi) ...
           - (tmp_x(pxind,:,:)-tmp_x(nxind,:,:))/(2*h) ...
           - (tmp_y(:,pyind,:)-tmp_y(:,nyind,:))/(2*h) ...
           - (tmp_z(:,:,pzind)-tmp_z(:,:,nzind))/(2*h) ;
       
    if (equation_type == 0)
        % Ginzburg-Landau
        dphidt = -M*dFdphi;        
    elseif (equation_type == 1)
        d2dFdphi = (dFdphi(pxind,:,:)+dFdphi(nxind,:,:) ...
                   +dFdphi(:,pyind,:)+dFdphi(:,nyind,:) ...
                   +dFdphi(:,:,pzind)+dFdphi(:,:,nzind)-6*dFdphi)/(h^2);
        dphidt = M*d2dFdphi;
    elseif (equation_type == 2)
        dphidt = -M*(dFdphi - mean(mean(mean(dFdphi))));
    else
        disp(sprintf('Unknown equation_type %d',equation_type));
        return
    end
    
    phi = phi + dt* dphidt;
    Gdata(iter) = max(max(max(abs(dphidt))));
    
    if(mod(iter,printfreq)==0)
        disp(sprintf('iter=%d F=%14.6e G=%14.6e',iter,F,Gdata(iter)));
    end
    
    if(mod(iter,plotfreq)==0)
        figure(1);
        clf
        p_surf = patch(isosurface(x,y,z,phi,0.5));
        isonormals(x,y,z,phi,p_surf)
        set(p_surf,'FaceColor','cyan','EdgeColor','none');
        daspect([1 1 1])
        view(3); axis tight
        camlight 
        lighting gouraud
        xlim([min(x) max(x)]); ylim([min(y) max(y)]); zlim([min(z) max(z)]);
        alpha(0.7);
        
        figure(2);
        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