Entry 3441
main.py
Submitted by name
on March 27, 2010 at 8:38 p.m.
Language: Python. Code size: 31.4 KB.
import numpy as np import sympy as sp from sympy.utilities import lambdify import matplotlib.pyplot as plt import pylab class DFT: def __init__(self, N = 128, dr = [0.085,0.085], sigma = 2.5, step_potential = True): """ Initializes a 2D DFT solver object. Inputs:- N is the number of points in the X and Y directions (128 default) dr is a tuple representing the grid spacing (set to 0.265 by default) step_potential is set to True to use the step potential model, set it to false to use archer's repulsive/attractive model. This will use dr = [0.265, 0.265] unless no_reset_dr is explicitly set to true Parameters that need to be set include:- mu The chemical potential T The temperature rho[:] The density field """ # Set the parameters self.sigma = sigma # length of step potential function self.N = N #N = 128 gives good results self.Dim = [N,N] self.step_potential = step_potential # comes before define_functional # i am forgetful and this will make sure I use the right dr value # for the potential model if not step_potential: print("Make sure you are using correct dr") self.dr = dr # changed to a list from a tuple because sometimes # it must be mutable self.define_functional() self.init_fields() self.set_wfs(dr) # Sets the weight functions with a given dr. #------------------------------------------------------------------------------ def define_functional(self): """ Defines the terms and weight functions which form the hard core FMT part as well as the mean field tail part. Also lambdifies these functions for use in the numerical solver.""" # defines the symbols that we use rho = sp.Symbol('rho') eta = sp.Symbol('eta') # eta is the packing fraction r = sp.Symbol('r') q = sp.Symbol('q') R = 0.5 # Just one component so hard core diameter fixed to 1.0 eta = sp.pi*rho*R**2.0 # 2D Packing fraction # Define the excess terms and their first and second derivitives self.Phi = rho*(sp.log(rho)-2.0-sp.log(1.0-eta)+(1.0-eta)**(-1.0)) self.Phi_l = lambdify('rho',self.Phi,"numpy") self.phi = sp.diff(self.Phi,rho) self.phi_l = lambdify('rho',self.phi,"numpy") if not self.step_potential: # Define the form of the mean field term: # Archers Double exponential version eps_a = eps_r = 1.0 # Could make these parameters. R_a = 1.0 R_r = 2.0 self.mfr = - eps_a/R_a**2.0*sp.exp(-r/R_a) \ + eps_r/R_r**2.0*sp.exp(-r/R_r) self.mfr_l = lambdify('r',self.mfr,"numpy") #------------------------------------------------------------------------------ def init_fields(self): """ sets up the fields""" self.rho = np.ones(self.Dim) self.c1 = np.zeros(self.Dim) self.mf = np.empty(self.Dim) self.mft = np.empty(self.Dim) self.vext = np.zeros(self.Dim) #------------------------------------------------------------------------------ def set_wfs(self,dr=None): """ (Re) Initilazes the weight functions """ if dr != None: if len(np.array(dr)) == 1: self.dr = (dr,dr) else: self.dr = dr # Want to make r position arrays. # Not very elegant way of making this position arrays. self.lr = np.empty([2,self.N]) for ax in range(2): self.lr[ax,:] = np.arange(0.0,self.dr[ax]*self.N,self.dr[ax]) self.lr[ax,self.N/2+1:] = np.flipud(self.lr[ax,1:self.N/2]) def r_func(i,j=0): x, y = self.lr[0,i], self.lr[1,j] return np.sqrt(x**2.0+y**2.0) self.r = np.fromfunction(r_func,self.Dim,dtype=int) if self.step_potential: # set up mean field with step potential # self.mf = self.dr[0]*self.dr[1]*np.fft.rfftn( # np.where(self.r<=self.sigma,1.0, # 0.0)) # set up mean field with step potential, this time including fudge # factor to account for aliasing artifacts caused by dr and sharp step # potential self.mf = np.where(self.r <= self.sigma, np.where( self.r == self.sigma, 0.5, 1.0), 0.0) * self.dr[0] * self.dr[1] print "scaling by,", (2.0*np.pi*self.sigma**2.0)/2.0/self.mf.sum() self.mf *= np.pi*2.5**2.0/self.mf.sum() self.mf = np.fft.rfftn(self.mf) else: # Added a limit just to make things a bit more formal! limit = 16.0 # Ideally could be longer! self.mf = self.dr[0]*self.dr[1]*np.fft.rfftn( np.where(self.r<limit,self.mfr_l(self.r)-self.mfr_l(limit), 0.0)) # Reset lr so that it is a bit more useful for ax in range(2): self.lr[ax,:] = np.arange(0.0,self.dr[ax]*self.N,self.dr[ax]) #------------------------------------------------------------------------------ def calc_c1(self,rho=None): """ Calculates all c1 including the log(rho) part. """ if(rho == None): rho = self.rho else: rho = rho.reshape(self.Dim) if (rho<0.0).any(): self.c1[:] = np.nan return self.c1.flatten() rhok = np.fft.rfftn(rho) self.mft = np.fft.irfftn(rhok*self.mf)/self.T self.c1[:] = np.where(rho>0.0, self.phi_l(rho) + self.mft + self.vext -self.mu ,0.0) return self.c1.flatten() #------------------------------------------------------------------------------ def calc_Omega(self,rho=None): """ Calculates the value of Omega """ if(rho == None): rho = self.rho else: rho = rho.reshape(self.Dim) # why? if (rho<0.0).any(): return 1000.0 rhok = np.fft.rfftn(rho) self.mft = np.fft.irfftn(rhok*self.mf)/self.T self.Omega = self.dr[0]*self.dr[1]*np.sum(np.where(rho>0.0, self.Phi_l(rho) + rho*(0.5*np.real(self.mft) + self.vext-self.mu),0.0)) return self.Omega #------------------------------------------------------------------------------ def set_mu(self,rhob,reset=False): """ Calculates the value of mu for a given rhob If reset = True then self.rho is set to rhob, else a temporary variable rho is set to rhob and self.rho is unchanged. In essence, setting reset = True causes the global self.rho value to be (re)set to rhob every time set_mu is run. NB: This function returns a value for mu but also sets the global value self.mu to the new mu value. It is a getter as well as a setter. Make sure vext[:] = 0.0 else it all goes wrong. """ if reset: rho = self.rho else: rho = np.empty(self.Dim) # why is this here paul? rho[:] = rhob self.mu = 0.0 self.mu = self.calc_c1(rho)[0] return self.mu #------------------------------------------------------------------------------ def Picard(self,rho=None,alpha=0.001,tol=1e-5): """ A simple Picard routine """ if(rho == None): rho = self.rho for outer in range(100): for inner in range(10): self.calc_c1(rho) rho -= alpha*self.c1 error = np.linalg.norm(self.c1.flatten(),ord=np.inf) print outer*10,error,self.calc_Omega(rho) if error < tol: break return rho #------------------------------------------------------------------------------ def jiggle(self,rho=None,a=-0.1,b=0.1): """ Mix up array rho a bit NB: I think we may be mixing it up too much """ if(rho == None): rho = self.rho rho += a+(b-a)*np.random.random(self.Dim) rho[:] = np.where(rho < 0.0, 0.000001, rho[:]) rho[:] = np.where(rho > 1.0,1.0, rho[:]) return rho #------------------------------------------------------------------------------ def Picard_cg(self,rho=None,alpha=0.04,tol=1e-9,i_count=100,visual = True, quiet = False): """ A conjugate gradient Picard routine Set visual to False to turn off plotting. """ # (re)set a flag in case we dont find a solution, we set this to 0 here # and in case of non-conversion we set it to 1 self.DNC=0 plt.ion() # Turns matplotlib to interactive mode if(rho == None): rho = self.rho k = 0 alpha_k = alpha beta_k = 0.0 gfk = self.calc_c1(rho) xk = rho.flatten() pk = -gfk gnorm = np.linalg.norm(gfk,ord=np.inf) if not quiet: print " k c1_norm Phi alpha beta" # k is not number of iterations, k is number of hundred iterations # under 100,000 would be nice, so lets set k<1000 while (gnorm > tol): if not quiet: print ' %5i, %10.6g %15.10g %10.6g %10.6g' \ % (k*i_count, gnorm, self.calc_Omega(xk), alpha_k, beta_k) k += 1 if visual == True: plt.clf() plt.imshow(xk.reshape(self.Dim), extent=(0,self.Dim[0]*self.dr[0],0,self.Dim[0]*self.dr[1])) plt.colorbar() plt.draw() for inner in xrange(i_count): deltak = np.dot(gfk,gfk) # alpha_k = alpha # Could change this to simple line search xk += alpha_k*pk gfkp1 = self.calc_c1(xk) yk = gfkp1 - gfk beta_k = max(0.0,np.dot(yk,gfkp1)/deltak) pk = -gfkp1 + beta_k * pk gfk = gfkp1 gnorm = np.linalg.norm(gfk,ord=np.inf) # Check if it has failed, if so reduce alpha and restart from beginning if( np.isnan(xk).any() == True or np.isnan(gnorm)): xk = rho.flatten() gfk = self.calc_c1(rho) pk = -gfk gnorm = np.linalg.norm(gfk,ord=np.inf) alpha_k /= 2.0 print "Adjusting alpha to %2.4g, beta was %2.4g" % \ (alpha_k, beta_k) # break if we hit our limit (25,000 iters) - if it hasnt got there # by then its taking too long anyway and probably is finding a # wrong solution if k >= 500: self.DNC = 1 # because we did not converge break # Non check for non-convergence ! rho[:] = xk.reshape(self.Dim) return rho # ----------------------------------------------------------------------------- def plot(self,array=None,save = False): """ Simply plots any array using matplotlib """ if array == None: array = self.rho # plt.close('all') plt.ioff() # Turns matplotlib off interactive mode. plt.imshow(array, extent=(0,self.Dim[0]*self.dr[0],0,self.Dim[1]*self.dr[1])); plt.colorbar() # plt.show() if save: filename = "./overview/T0.1/" + "rho_" + str(np.average(self.rho))[0:5] + ".eps" plt.savefig(filename) #------------------------------------------------------------------------------ def add_ext_pot(self,pot=None,params=(0.0,1.0)): """ Add external potentials. At moment only adds a particle of the type specified in the model to the corner of the region. """ if pot == "Zero": self.vext[:] = 0.0 if pot == "Particle": self.vext[:] += self.mfr_l(self.r) self.vext[:] += np.where(self.r<=1.0,np.inf,0.0) # Apply ideal gas approximation to rho. self.rho[:] *= np.exp(-self.vext) #------------------------------------------------------------------------------ def fastrun(self, T = 0.2, rhob = 0.2, quiet = True): """ A simple running of the program that completes fairly fast """ self.T = T self.set_mu(rhob=rhob,reset=True) self.jiggle() self.Picard_cg(alpha=0.04,tol=1e-6,visual=True,quiet=quiet) #------------------------------------------------------------------------------ def save_rho(self,file_name): """Simple function to save rho """ res = np.savetxt(file_name,self.rho) return res #------------------------------------------------------------------------------ def load_rho(self,file_name): """Simple function to save rho """ self.rho = np.loadtxt(file_name) return self.rho #----------------------------------------------------------------------------- def test(self,rhob,T): """ A simple use of the program, identical to the example use, but it takes bulk density (rhob) and temperature (T) arguments. """ # make input variables into objects self.T = T # i think reset = True must be set, not sure why yet self.set_mu(rhob=rhob,reset=True) # jiggle() obviously randomises the density grid self.jiggle() # temporarily set alpha to a lower value here self.Picard_cg(alpha=0.005,tol=1e-6) # picard routine plots graphs as it goes so im assuming this just plots # the final one self.plot() #------------------------------------------------------------------------------ def solve(self, rhob ,T ,alpha=0.02 ,tol=1e-8, reset = False): """ A simple solve function intended for use by scripts and batch runs If reset = True then density grid is reset to bulk density after every run. Else, it just keeps the old one around. """ # make input variables into objects self.T = T self.set_mu(rhob=rhob,reset=reset) # No need to jiggle density grid if it hasn't been reset, as it is # already seeded. This may cause problems if the program hasn't been # run yet and the density grid is not already defined. Although this # may be ok if the set_mu function compensates for that. Either way if # reset = False it would be a good idea to solve a "seeder" state point # first in order to provide a base to build on. if reset: self.jiggle(a=0.01,b=-0.01) # temporarily set alpha to a lower value here self.Picard_cg(alpha=alpha,tol=tol) # i do not think we need to plot anything # self.plot() #------------------------------------------------------------------------------ def get_omega_density(self): """ Get grand potential (omega) and divide by area of the grid to obtain an energy density (pressure). This can be useful for synchronising fluid modulations to the grid size in order to avoid forcing the period to be Lx/n in the x direction and Ly/m in the y direction (where n and m are integers). """ omega = self.calc_Omega() # calculate lengths of grid l_x = self.Dim[0] * self.dr[0] l_y = self.Dim[1] * self.dr[1] area = l_x * l_y omega_density = omega / area print("inside get_omega_density function:") print("omega = " + str(omega) + ", dr = " + str(self.dr) + ", omega_density = " + str(omega_density)) return omega_density #------------------------------------------------------------------------------ def plot_dr(self, reset = True): """ Iterate through grid spacing values from start_dr to end_r, call get_omega_density for each and plot values in a list, then return numpy array of results. increment_dr is the step size taken each time This gets a lot harder with two independent values for dr[0] and dr[1] We will use identical values for each as this seems to usually be a minima anyway, and even if it isn't it's close enough for our purposes We will use T=0.2 and rhob = 0.2 just because these are fairly standard values NB: We should probably not reset the density grid each time otherwise calculation will get difficult and it will take a long time as well as possibly finding different minima each time. This would cause problems. """ # create two dimensional numpy array to put results in self.omega_array = [[], []] # this seems to be a stable state point self.T = 0.2 self.rhob = 0.35 # this is here in order to initialise a sample density grid before we # try and solve for a uniform one and get the wrong answer for i in range(0, 120):# cos end point is never returned by range() self.dr[0] = 0.08 + (float(i)/10000) self.dr[1] = self.dr[0] # keep dr_x and dr_y identical self.set_wfs(self.dr) print("i = " + str(i)) print("self.dr = " + str(self.dr)) # solve for the given state point self.solve(T = self.T, rhob = self.rhob, reset = reset, tol=1e-10) omega_dens = self.get_omega_density() print("omega_dens = " + str(omega_dens)) # put values into "array", omega_dens[0] is dr, omega_dens[1] is mu # density self.omega_array[0].append(self.dr[0]) self.omega_array[1].append(omega_dens) # a pretty(ish) printout # print("dr omega_dens") # for i in range(0, 121): # print("%.3f %.6g" % (self.omega_array[0][i], # self.omega_array[1][i])) # let's export a numpy array omega_array = np.array(self.omega_array) return omega_array #------------------------------------------------------------------------------ def get_overview(self, T, rhob_low = 0.05, rhob_high = 0.8, rhob_step = 0.02, reset = False, save = False, phase = None, rho = False): """Iterate through a combination of every rhob and T value given, and output a list of rhob and corresponding other values for each T, in one big list. Arguments for rhob are self-explanatory, T is the temperature. If reset is true then function will solve from a random grid for each point. Else, it will use the grid from last time. Supplying a negative rhob will cause the function to automatically plot backwards instead of forwards. NB: If running backwards, rhob_start and rhob_end should be supplied to the function reversed (start high, get smaller with negative rhob_step, end low). This is the logical way to do it (I think). A density array can be provided as a starting point. This is given with the phase variable which is a string, either None (for fluid), "stripes", "blobs", or "clusters". A bit of hacked up string building should automagically load the correct text file array for the given T and phase. If a density array (sys.rho) is supplied as an argument to rho, phase will be ignored and the supplied density grid used instead. """ ## Initialise variables iters = 0 # are we running backwards or forwards? if rhob_step < 0: backwards = True rhob_start = rhob_high rhob_end = rhob_low else: backwards = False rhob_start = rhob_low rhob_end = rhob_high # Use while loops instead of for loops because although its less # elegant. It is also simpler and I can't be arsed mucking around # converting floats to integers for the ranges. # Set up array # [[rhob], [rho_ave], [chem_pot], [pressure]] pressure_array = [[], [], [], []] self.T = T # set up initial rhob rhob = rhob_start # do a setup solve IF rho is not set (if rho is set we use old sys.rho) if rho == False: if phase == None: self.solve(rhob = rhob, T = T, reset = True, tol = 1e-8) elif phase != None: # some filename loading logic path = "../saved_rho/machine_readable/" filename = "T" + str(T) + "_" + phase self.rho = np.loadtxt(path + filename) # I don't like this code duplication, it seems fundamentally wrong, but # it's the quick and dirty way to make this work. I probably should # have used for loops, or put all this stuff in a function somewhere. # Oh well, it works. if backwards: while rhob > rhob_low: iters += 1 print "This is iteration number " + str(iters) print "T = " + str(T) print "rhob = " + str(rhob) self.rhob = rhob # Call solver function self.solve(rhob = rhob, T = T, reset = reset, tol = 1e-9) # a little jiggle to knock us out of false minima? probably not # necessary and may even hurt things # self.jiggle(a=0.001,b=-0.001) # append rhob pressure_array[0].append(rhob) # append rho_ave pressure_array[1].append(np.average(self.rho)) # append chemical potential pressure_array[2].append(self.mu) # append pressure pressure_array[3].append(self.get_omega_density()) # to save a picture: if save: plt.clf() plt.imshow(self.rho, extent=(0,self.Dim[0]*self.dr[0],0,self.Dim[1]*self.dr[1])); plt.colorbar() plt.draw() filename = "./overview/T" + str(T) + "/" + "rho_" +\ str(np.average(self.rho))[0:5] + "_backward.eps" plt.savefig(filename) # Remember to increment rhob rhob += rhob_step else: while rhob < rhob_high: iters += 1 print "This is iteration number " + str(iters) print "T = " + str(T) print "rhob = " + str(rhob) self.rhob = rhob # Call solver function self.solve(rhob = rhob, T = T, reset = reset, tol = 1e-9) # a little jiggle to knock us out of false minima? probably not # necessary and may even hurt things # self.jiggle(a=0.001,b=-0.001) # append rhob pressure_array[0].append(rhob) # append rho_ave pressure_array[1].append(np.average(self.rho)) # append chemical potential pressure_array[2].append(self.mu) # append pressure pressure_array[3].append(self.get_omega_density()) # to save a picture: if save: plt.clf() plt.imshow(self.rho, extent=(0,self.Dim[0]*self.dr[0],0,self.Dim[1]*self.dr[1])); plt.colorbar() plt.draw() filename = "./overview/T" + str(T) + "/" + "rho_" +\ str(np.average(self.rho))[0:5] + "_forward.eps" plt.savefig(filename) # Remember to increment rhob rhob += rhob_step # let's make it a numpy array pressure_array = np.array(pressure_array) # let's save it to be sure if self.step_potential == True: filename = "step_" + "T_" + str(self.T) + "rhob_" +\ str(rhob_start) + "-" + str(rhob_end) + ".txt" else: filename = "arch_" + "T_" + str(self.T) + ".txt" np.savetxt(filename, pressure_array) return pressure_array #------------------------------------------------------------------------------ def traverse_T(self, T_low, T_high, rhob, T_step = 0.02, reset = False, save = False, rho = False): """ This is an analog of the get_overview function but for set rhob and varying T. Iterate through a combination of every T value for a given rhob, and output a list of T and corresponding other other values in one big list. Arguments for T are self-explanatory, rhob is bulk density. If reset is true then function will solve from a random grid for each point. Else, it will use the grid from last time. Supplying a negative rhob will cause the function to automatically plot downwards instead of upwards (downwards = decreasing T). NB: If running downwards, T_start and T_end should be supplied to the function reversed (start high, get smaller with negative T_step, end low). This is the logical way to do it (I think). Setting rho to true will use the old density grid (sys.rho) from the previous solve """ ## Initialise variables iters = 0 # are we running upward or downward? if T_step < 0: downward = True T_start = T_high T_end = T_low else: downward = False T_start = T_low T_end = T_high # Use while loops instead of for loops because although its less # elegant. It is also simpler and I can't be arsed mucking around # converting floats to integers for the ranges. # Set up array # [[T], [rho_ave], [chem_pot], [pressure], [rhob]] pressure_array = [[], [], [], [], []] # set up initial T T = T_start # do a setup solve IF rho is not set (if rho is set we use old sys.rho) if rho == False: self.solve(rhob = rhob, T = T, reset = True, tol = 1e-8) # I don't like this code duplication, it seems fundamentally wrong, but # it's the quick and dirty way to make this work. I probably should # have used for loops, or put all this stuff in a function somewhere. # Oh well, it works. if downward: while T > T_low: iters += 1 print "This is iteration number " + str(iters) print "rhob = " + str(rhob) print "T = " + str(T) self.T = T # Call solver function self.solve(rhob = rhob, T = T, reset = reset, tol = 1e-9) # a little jiggle to knock us out of false minima? probably not # necessary and may even hurt things # self.jiggle(a=0.001,b=-0.001) # append T pressure_array[0].append(T) # append rho_ave pressure_array[1].append(np.average(self.rho)) # append chemical potential pressure_array[2].append(self.mu) # append pressure pressure_array[3].append(self.get_omega_density()) # append rhob (not strictly necessary as we can solve from mu # but nice to have anyway pressure_array[4].append(rhob) # to save a picture: if save: plt.clf() plt.imshow(self.rho, extent=(0,self.Dim[0]*self.dr[0],0,self.Dim[1]*self.dr[1])); plt.colorbar() plt.draw() filename = "./overview/rhob" + str(rhob) + "/" + "T_" +\ str(self.T) + ".eps" plt.savefig(filename) # Remember to increment T T += T_step else: while T < T_high: iters += 1 print "This is iteration number " + str(iters) print "rhob = " + str(rhob) print "T = " + str(T) self.T = T # Call solver function self.solve(rhob = rhob, T = T, reset = reset, tol = 1e-9) # a little jiggle to knock us out of false minima? probably not # necessary and may even hurt things # self.jiggle(a=0.001,b=-0.001) # append T pressure_array[0].append(T) # append rho_ave pressure_array[1].append(np.average(self.rho)) # append chemical potential pressure_array[2].append(self.mu) # append pressure pressure_array[3].append(self.get_omega_density()) # append rhob (not strictly necessary as we can solve from mu # but nice to have anyway pressure_array[4].append(rhob) # to save a picture: if save: plt.clf() plt.imshow(self.rho, extent=(0,self.Dim[0]*self.dr[0],0,self.Dim[1]*self.dr[1])); plt.colorbar() plt.draw() if downward: filename = "./overview/rhob" + str(rhob) + "/d" + "T_" +\ str(self.T) + ".eps" else: filename = "./overview/rhob" + str(rhob) + "/u" + "T_" +\ str(self.T) + ".eps" plt.savefig(filename) # Remember to increment T T += T_step # let's make it a numpy array pressure_array = np.array(pressure_array) # let's save it to be sure if self.step_potential == True: filename = "step_" + "rhob_" + str(rhob) + "T_" +\ str(T_start) + "-" + str(T_end) + ".txt" else: filename = "arch_" + "rhob_" + str(rhob) + "T_" +\ str(T_start) + "-" + str(T_end) + ".txt" np.savetxt(filename, pressure_array) return pressure_array #------------------------------------------------------------------------------ if __name__ == "__main__": # an example of how to instantiate the program example = DFT(128,[0.25,.25]) example.test(0.2,0.2)
This snippet took 0.11 seconds to highlight.
Back to the Entry List or Home.