% TDGL_A_2D.m % This MATLAB program simulates the dynamics of the % Landau phi^4 free energy model in 2-D using % Ginzburg-Landau TDGL model A dynamics. % % K. Beers. MIT ChE. 1/16/04. v 7/12/05 function iflag_main = TDGL_A_2D(); iflag_main = 0; disp('Running TDGL_A_2D : '); disp('Sample default values are in () in input prompts.'); % Ask the user to specify the length of the % system. disp(' '); disp('Enter system parameters : '); Grid.L = input('----- Enter side length of simulation domain (100) : '); Grid.V = Grid.L^2; % the total 2-D "volume" % Next, ask the user to specify the number of % grid points in each direction of the 2D grid. Grid.N_side = input('----- Enter # of grid points in each direction (50) : '); Grid.N_tot = Grid.N_side^2; Grid.x = linspace(0,Grid.L,Grid.N_side); Grid.y = linspace(0,Grid.L,Grid.N_side); Grid.d = 2; % dimension of the lattice % compute the grid spacing and the 2-D volume per cell Grid.h = Grid.L/Grid.N_side; % we use periodic BC Grid.v_cell = Grid.h^2; % To parameterize the Landau model, we use the Ising lattice % model in 2-D. We define energy scale such that % kb*Tc = J*z = 1. Sim.k_b = 1; % Boltzmann constant Sim.z = 2*Grid.d; % # of nearest neighbors per site Sim.J = 1/Sim.z; % set coupling constant such that J*z = 1 Sim.Tc_MF = Sim.J*Sim.z/Sim.k_b; % mean-field critical temperature Sim.Tc_RG = Sim.J/Sim.k_b/0.44069; % renormalization-group Tc % We now ask the user to input the temperature disp(' '); disp('Next, the temperature is entered, in comparison to : '); disp([' mean-field Tc = ', num2str(Sim.Tc_MF)]); disp([' Renormalization group Tc = ', num2str(Sim.Tc_RG)]); Sim.T = input('----- Enter simulation T (0.35) : '); % enter field strength in units of kb*Tc disp(' '); disp('H is the external field strength, spin moment = 1.'); var1 = input('----- Enter H in units of kb*Tc (0) : '); Sim.H = var1*Sim.k_b*Sim.Tc_MF; % We now compute the parameters of the Landau % phenomenological free energy for this Ising % lattice. This is obtained by expanding the % free energy around phi = 0, and retaining up % to the fourth power terms (leading error is % then 6th power due to symmetry). Sim.r0 = Sim.k_b/Grid.v_cell; Sim.r = Sim.r0*(Sim.T-Sim.Tc_MF); Sim.u = Sim.k_b*Sim.T/12/Grid.v_cell; Sim.c = Sim.J*Grid.h^2/Grid.v_cell; Sim.w = Sim.H / Grid.v_cell; % Ask the user to input the transport coefficient disp(' '); disp('Gamma is the transport coefficient.'); Sim.gamma = input('----- Enter transport coefficient Gamma (1) : '); % We now ask the user to input the time step of the simulation, % after reporting an approximate upper limit time step of an % explicit integrator for a 2-D diffusion problem. Since % an implicit algorithm is used, this time step limitation % is not enforced, however. Sim.dt_max_explicit = (Grid.h^2)/(Sim.gamma*Sim.c); disp(' '); disp('Next, you are asked to set time step, in comparison to : '); disp([' Explicit integrator critical dt (upper bound) = ', ... num2str(Sim.dt_max_explicit)]); Sim.dt = input('----- Enter time step for implicit integrator (0.5) : '); Sim.t_end = input('----- Enter the end time for the simulation (2000) : '); Sim.num_steps = 1 + floor(Sim.t_end / Sim.dt); Sim.plot_freq = input('----- Enter # of time steps between plots (10) : '); % Generate std of random noise Sim.rand_noise_std = sqrt(2*Sim.k_b*Sim.T*Sim.gamma/Grid.v_cell*Sim.dt); % We now allocate a vector to store the field values at each % grid point. We initialize this vector to store the initial % values of all zero. phi = zeros(Grid.N_tot,1); % We next specify a numerical stencil operator that approximates % the negative 2-D Laplacian on a grid with unit spacing. Here, we use the % traditional finite difference formula; however, a better choice would % be to include corner sites to obtain a more isotropic stencil operator % to lessen the effect of the lattice orientation upon the results. Grid.A = spalloc(Grid.N_tot,Grid.N_tot,5*Grid.N_tot); for ix=1:Grid.N_side for iy=1:Grid.N_side % center site l = get_label_2D(ix,iy,Grid.N_side); Grid.A(l,l) = 4; % -x neighbor jx = ix-1; jy = iy; if(jx < 1) % periodic BC jx = jx + Grid.N_side; end k = get_label_2D(jx,jy,Grid.N_side); Grid.A(l,k) = -1; % +x neighbor jx = ix+1; jy = iy; if(jx > Grid.N_side) jx = jx - Grid.N_side; end k = get_label_2D(jx,jy,Grid.N_side); Grid.A(l,k) = -1; % -y neighbor jx = ix; jy = iy - 1; if(jy < 1) jy = jy + Grid.N_side; end k = get_label_2D(jx,jy,Grid.N_side); Grid.A(l,k) = -1; % +y neighbor jx = ix; jy = iy + 1; if(jy > Grid.N_side) jy = jy - Grid.N_side; end k = get_label_2D(jx,jy,Grid.N_side); Grid.A(l,k) = -1; end end % We next compute the initial chemical potential field mu = calc_chem_pot(phi,Grid,Sim); % Finally, we open a figure that shows the order parameter % field dynamics. fig_order = figure; % and add figure for plot of average order parameter vs. time fig_order_mean = figure; % We now begin time simulation. time = 0; iframe = 0; for itime = 1:Sim.num_steps time = time + Sim.dt; % generate random noise for each grid point rand_noise = Sim.rand_noise_std*randn(Grid.N_tot,1); % get new field values using explicit Euler method phi = phi - Sim.gamma*Sim.dt*mu + rand_noise; % compute average order parameter phi_mean = mean(phi); % compute new chemical potential field mu = calc_chem_pot(phi,Grid,Sim); % periodically print out the results of the simulation to % the screen, for use in making a movie. if(~mod(itime,Sim.plot_freq)) iframe = iframe + 1; figure(fig_order); % extract results into 2-D array for plotting Phi_plot = zeros(Grid.N_side,Grid.N_side); for ix=1:Grid.N_side for iy=1:Grid.N_side l = get_label_2D(ix,iy,Grid.N_side); Phi_plot(ix,iy) = phi(l); end end phi_max = norm(phi,inf); contourf(Grid.x,Grid.y,Phi_plot); title_phrase = 'TDGL-A \phi^4 model: '; title_phrase = [title_phrase, 'T = ', num2str(Sim.T)]; title_phrase = [title_phrase, ', H = ', num2str(Sim.H)]; title_phrase = [title_phrase, ', time = ', num2str(time)]; xlabel('x'); ylabel('y'); title(title_phrase); colorbar; Movie_order(iframe) = getframe; % add plot of average order parameter vs. time figure(fig_order_mean); hold on; plot(time,phi_mean,'o'); xlabel('time'); ylabel('< \phi >'); title(title_phrase); end end % show movie and convert to avi repeat_movie = 1; fig_movie = figure; while(repeat_movie) figure(fig_movie); movie(Movie_order,0); repeat_movie = input('Repeat movie (0=n, 1=y) : '); end iflag_main = 1; return; % -------------------------------------------- % This routine returns an index label for each % point in a 2-D grid, based on the same % numbering system used in meshgrid and the % plotting of results. The x-axis goes from % left to right, and the y-axis from top to % bottom, in the same ordering as the elements % of a matrix. function label = get_label_2D(ix,iy,N_side); label = (ix-1)*N_side + iy; return; % --------------------------------------------- % This routine returns the real-space coordinates % of a grid point in the 2-D lattice, using the % same numbering system as meshgrid. (x_O,y_O) % are the coordinates of the origin, located at the % upper left corner of the domain. function [x_c,y_c] = get_coord_2D(ix,iy,h,x_O,y_O); x_c = x_0 + h*(ix-1); y_c = y_O - h*(iy-1); return; % ------------------------------------------------------ % This routine computes the chemical potential field % for the Landau phi^4 model on a 2-D lattice grid. function mu = calc_chem_pot(phi,Grid,Sim); mu = Sim.c/(Grid.h^2)*Grid.A*phi + Sim.r*phi + ... 4*Sim.u*phi.^3 - Sim.w*ones(size(phi)); return;