%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%File: pf3d_single_iso.m
%
%Purpose: 3-dimensional Phase Field model 
%         solving Ginzburg-Landau (or Cahn-Hillard) equation
%         by finite difference method
%        
% Wei Cai (caiwei@stanford.edu)
% Seunghwa Ryu (shryu@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'))
    model_param = [0.4 1.0 1.0 1.0]; 
end

h = model_param(1);      % grid size
M = model_param(2);      % kinetic coefficient
U = model_param(3);      % energy parameters
vareps = model_param(4); % energy parameters

if(~exist('Niter'))
    Niter = 5000; % 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 = rand(N_grid);
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);

    F = (U * sum(sum(sum( phi.^2.*(1-phi).^2 ))) ...
        + vareps * sum(sum(sum( dphidx.^2 + dphidy.^2 + dphidz.^2)))) * (h^3);
    Fdata(iter)=F;
   
    dFdphi = 2*U*phi.*(1-phi).*(1-2*phi) - 2*vareps * d2phi;
 
    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