Révision 88 Ising/MP/Ising2D-MP.py
| Ising2D-MP.py (revision 88) | ||
|---|---|---|
| 33 | 33 |
# Random access coordonate |
| 34 | 34 |
X,Y=numpy.random.randint(SizeX),numpy.random.randint(SizeY) |
| 35 | 35 |
|
| 36 |
DeltaE=J*sigma[X,Y]*(2*(sigma[X,(Y+1)%SizeY]+
|
|
| 36 |
DeltaE=sigma[X,Y]*(2*J*(sigma[X,(Y+1)%SizeY]+
|
|
| 37 | 37 |
sigma[X,(Y-1)%SizeY]+ |
| 38 | 38 |
sigma[(X-1)%SizeX,Y]+ |
| 39 | 39 |
sigma[(X+1)%SizeX,Y])+B) |
| ... | ... | |
| 46 | 46 |
def Magnetization(sigma,M): |
| 47 | 47 |
return(numpy.sum(sigma)/(sigma.shape[0]*sigma.shape[1]*1.0)) |
| 48 | 48 |
|
| 49 |
def Energy(sigma,J): |
|
| 49 |
def Energy(sigma,J,B):
|
|
| 50 | 50 |
# Copier et caster |
| 51 | 51 |
E=numpy.copy(sigma).astype(numpy.float32) |
| 52 | 52 |
|
| 53 | 53 |
# Appel par slice |
| 54 |
E[1:-1,1:-1]=-J*E[1:-1,1:-1]*(E[:-2,1:-1]+E[2:,1:-1]+
|
|
| 55 |
E[1:-1,:-2]+E[1:-1,2:])
|
|
| 54 |
E[1:-1,1:-1]=E[1:-1,1:-1]*(2.*J*(E[:-2,1:-1]+E[2:,1:-1]+
|
|
| 55 |
E[1:-1,:-2]+E[1:-1,2:])+B)
|
|
| 56 | 56 |
|
| 57 | 57 |
# Bien nettoyer la peripherie |
| 58 | 58 |
E[:,0]=0 |
| ... | ... | |
| 64 | 64 |
|
| 65 | 65 |
return(Energy/(E.shape[0]*E.shape[1]*1.0)) |
| 66 | 66 |
|
| 67 |
def CriticalT(T,E): |
|
| 68 |
|
|
| 69 |
Epoly=numpy.poly1d(numpy.polyfit(T,E,T.size/3)) |
|
| 70 |
dEpoly=numpy.diff(Epoly(T)) |
|
| 71 |
dEpoly=numpy.insert(dEpoly,0,0) |
|
| 72 |
return(T[numpy.argmin(dEpoly)]) |
|
| 73 |
|
|
| 67 | 74 |
def DisplayCurves(T,E,M,J,B): |
| 68 | 75 |
|
| 69 | 76 |
plt.xlabel("Temperature")
|
| ... | ... | |
| 99 | 106 |
try: |
| 100 | 107 |
opts, args = getopt.getopt(sys.argv[1:],"hcj:b:z:i:s:e:p:u:",["coupling=","magneticfield=","size=","iterations=","tempstart=","tempend=","tempstep=","units"]) |
| 101 | 108 |
except getopt.GetoptError: |
| 102 |
print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -u <Unit of process> -c (Print Curves)' % sys.argv[0] |
|
| 109 |
print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Square Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -u <Unit of process> -c (Print Curves)' % sys.argv[0]
|
|
| 103 | 110 |
sys.exit(2) |
| 104 | 111 |
|
| 105 | 112 |
for opt, arg in opts: |
| 106 | 113 |
if opt == '-h': |
| 107 |
print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -c (Print Curves)' % sys.argv[0] |
|
| 114 |
print '%s -j <Coupling Factor> -b <Magnetic Field> -z <Size of Square Lattice> -i <Iterations> -s <Minimum Temperature> -e <Maximum Temperature> -p <steP Temperature> -c (Print Curves)' % sys.argv[0]
|
|
| 108 | 115 |
sys.exit() |
| 109 | 116 |
elif opt == '-c': |
| 110 | 117 |
Curves=True |
| ... | ... | |
| 142 | 149 |
|
| 143 | 150 |
Trange=numpy.arange(Tmin,Tmax+Tstep,Tstep) |
| 144 | 151 |
|
| 145 |
E=[] |
|
| 146 |
M=[] |
|
| 147 |
|
|
| 148 | 152 |
def MetropolisStrip(T): |
| 149 | 153 |
# Indispensable d'utiliser copy : [:] ne fonctionne pas avec numpy ! |
| 150 | 154 |
sigma=numpy.copy(sigmaIn) |
| 151 | 155 |
duration=Metropolis(sigma,J,B,T,Iterations) |
| 152 |
ImageOutput(sigma,"Ising2D_Serial_%i_%1.1f_Final" % (Size,T))
|
|
| 156 |
ImageOutput(sigma,"Ising2D_MP_%i_%1.1f_Final" % (Size,T))
|
|
| 153 | 157 |
print "CPU Time : %f" % (duration) |
| 154 |
E,M=Energy(sigma,J),Magnetization(sigma,B) |
|
| 158 |
indice=numpy.where(Trange==T)[0][0] |
|
| 159 |
E,M=Energy(sigma,J,B),Magnetization(sigma,B) |
|
| 155 | 160 |
print "Total Energy at Temperature %f : %f" % (T,E) |
| 156 | 161 |
print "Total Magnetization at Temperature %f : %f" % (T,M) |
| 157 |
return(T,E,M)
|
|
| 162 |
return([T,E,M])
|
|
| 158 | 163 |
|
| 159 | 164 |
pool=Pool(processes=Procs) |
| 160 | 165 |
print "Start on %i processors..." % Procs |
| ... | ... | |
| 162 | 167 |
# MAP: distribution of usecases T to be applied to MetropolisStrip |
| 163 | 168 |
Results=pool.map(MetropolisStrip,Trange) |
| 164 | 169 |
|
| 170 |
E=numpy.array(Results)[:,1] |
|
| 171 |
M=numpy.array(Results)[:,2] |
|
| 172 |
|
|
| 165 | 173 |
if Curves: |
| 166 | 174 |
DisplayCurves(Trange,E,M,J,B) |
| 167 | 175 |
|
| 168 | 176 |
# Save output |
| 169 | 177 |
numpy.savez("Ising2D_MP_%i_%.8i" % (Size,Iterations),(Trange,E,M))
|
| 170 |
|
|
| 178 |
|
|
| 179 |
# Estimate Critical temperature |
|
| 180 |
print "The critical temperature on %ix%i lattice with J=%f, B=%f is %f " % (Size,Size,J,B,CriticalT(Trange,E)) |
|
| 181 |
|
|
Formats disponibles : Unified diff