Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
- #include <sys/stat.h>
- // =======================================================================
- // Global variables and helper functions you can use
- // -----------------------------------------------------------------------
- signed char spins[128][128]; //!< The spin array: either -1 (down) or 1 (up)
- int extent; //!< The extent (size) of the lattice
- float beta; //!< \beta
- int num_iter; //!< number of iterations
- double rnd(); //!< returns a random number in [0,1]
- int rnd_spin(); //!< returns a random spin, either -1 or 1
- //! Used by the test routine, checks if you set the next neighbours correctly
- void check_nn(int, int, int, int, int, int, int) ;
- // =======================================================================
- // Declarations of all functions you need to edit
- // -----------------------------------------------------------------------
- void initialize(int); //!< initialize spins
- void sweep(); //!< run an update sweep
- void measure(); //!< measure magnetisation
- int calculate_spin(int, int, int, int, int); //!< return a new spin
- // =======================================================================
- // A routine to initalize the spin array
- // -----------------------------------------------------------------------
- void initialize(int start_cond) {
- // Loop through the spin array
- int x,y;
- for(x=0; x<extent; x++)
- for(y=0; y<extent; y++) {
- // /////////////////////////////////////////////////////////////////////////
- /* Task 1: To warm up, something easy first: Depending on the starting
- * condition (variable "start_cond"), the spins should be
- * initialized:
- *
- * 0: random spins, 1: all spins up (1), 2: all spins down (-1)
- *
- * A loop over the lattice is alreaddy there, you only need
- * to check the starting condition (e.g. with an if-statement)
- * and set spin[x][y] accordingly. Use rnd_spin() to get random
- * spins....
- */
- // your code here ...
- if (start_cond == 0) {
- spins[x][y] = rnd_spin();
- }
- else if (start_cond == 1) {
- spins[x][y] = 1;
- }
- else if (start_cond == 2) {
- spins[x][y] = -1;
- }
- // /////////////////////////////////////////////////////////////////////////
- }
- }
- // =======================================================================
- // A routine to update the spin array
- // -----------------------------------------------------------------------
- void sweep() {
- // variables
- int cur; // the current spin
- int lns, rns, uns, dns; // left, right, up and down neighbour to the spin
- int x,y,eo; // loop variables
- for (eo=0;eo<=1;eo++) // even/odd splitting
- for (x=0;x<extent;x++) // loop x
- for (y=0;y<extent;y++) { // loop y
- // Do nothing if the coordinate (x,y) we are currently
- // on does not match the even/odd variable
- if (((x+y)%2)!=eo) continue;
- // get the current spin
- cur = spins[x][y];
- // /////////////////////////////////////////////////////////////////////////
- /* Task 2: This is inside a loop over the whole lattice. The variable "cur"
- * now contains the spin at x,y. You must fill the variables
- * "lns", "rns", "uns", "dns" with the correct neighbours:
- * left: x-1, right: x+1, up: y-1, down: y+1
- *
- * lns = ....
- * rns = ....
- *
- * Remember that there are boundary conditions, so take care
- * for x,y=0 and x,y=(extent-1)
- */
- // your code here
- if (x > 0 && x < (extent-1) && y > 0 && y < (extent-1)) {
- lns = spins[x-1][y];
- rns = spins[x+1][y];
- uns = spins[x][y-1];
- dns = spins[x][y+1];
- }
- else if (y == 0 && x > 0 && x < (extent-1)) {
- lns = spins[x-1][y];
- rns = spins[x+1][y];
- dns = spins[x][y+1];
- }
- else if (y == (extent-1) && x > 0 && x < (extent-1)) {
- lns = spins[x-1][y];
- rns = spins[x+1][y];
- uns = spins[x][y-1];
- }
- else if (x == 0 && y > 0 && y < (extent-1)) {
- rns = spins[x+1][y];
- uns = spins[x][y-1];
- dns = spins[x][y+1];
- }
- else if (x == (extent-1) && y > 0 && y < (extent-1)) {
- lns = spins[x-1][y];
- uns = spins[x][y-1];
- dns = spins[x][y+1];
- }
- else if (x == 0 && y == 0) {
- rns = spins[x+1][y];
- dns = spins[x][y+1];
- }
- else if (x == 0 && y == (extent-1)) {
- rns = spins[x+1][y];
- uns = spins[x][y-1];
- }
- else if (x == (extent-1) && y == 0) {
- lns = spins[x-1][y];
- dns = spins[x][y+1];
- }
- else if (x == (extent-1) && y == (extent-1)) {
- lns = spins[x-1][y];
- uns = spins[x][y-1];
- }
- // /////////////////////////////////////////////////////////////////////////
- // if you compile the test, this checks if you set the neighbours right
- // if you compile the normal programm, this does nothing
- check_nn(x,y,eo,lns,rns,uns,dns);
- // update the current spin (you will have to complete Task 3):
- spins[x][y] = calculate_spin(cur, lns, rns, uns, dns);
- }
- }
- // =======================================================================
- // Calculate a new spin from its neighbours and beta
- // -----------------------------------------------------------------------
- int calculate_spin(int cur, int lns, int rns, int uns, int dns) {
- int new_spin;
- // ////////////////////////////////////////////////////////////////////////
- /* Task 3 : Now some physics:
- * Here you have to calculate a new spin. "cur" holds the spin, "lns",
- * "rns", "uns" and "dns" its neighbours (o.k., you probably know that
- * because you just wrote that part ....)
- *
- * As you (should) know, in the Ising model the spin flips if the flipped
- * position is energetically favored or equal. If it energetically
- * disfavored, there is still a (temperature-dependent) chance
- * the spin will flip:
- *
- * P_flip = exp(-beta * (energy difference))
- *
- * You can (and should) use "rnd()" to generate a random number
- * between 0 and 1 to implement this.
- *
- */
- // your code here ...
- // /////////////////////////////////////////////////////////////////////////
- return new_spin;
- }
- // =======================================================================
- // Measures and outputs energy and magnetisation to file
- // (it's already complete, you do not need to edit it)
- // -----------------------------------------------------------------------
- void measure(FILE* output, int iteration) {
- int x,y;
- double ene=0; // Store enery and magnetisation
- double mag=0;
- // Energy via the Hamiltonian - no boundary-crossing terms,
- // but on average they don't matter
- for(y=0;y<extent-1;y++)
- for(x=0;x<extent-1;x++)
- ene += 2*spins[x][y] * ( spins[x][y+1] + spins[x+1][y] );
- // Magnetisation as sum of spin directions
- for(y=0;y<extent;y++)
- for(x=0;x<extent;x++)
- mag += spins[x][y];
- // Both have to be divided, so they are "per site" values
- ene /= (extent-1)*(extent-1);
- mag /= extent*extent;
- // Write the result to the output file
- fprintf(output, "%6d\t%e\t%e\n", iteration, ene, fabs(mag));
- }
- /* ****************** END ************************************************
- *
- * Below here, the main routine that controls the programm and the
- * helper functions are implemented. You do not need to change anything
- * here, you only need to read these if you are interested.
- *
- * ********************************************************************** */
- #ifndef _TEST_RUN //Used for test-ising, don't modify these "#ifndef" things
- //! The main routine controlling program flow
- int main(int argc, char** argv) {
- // Check the number of command line arguments
- if ((argc!=4)&&(argc!=5)) {
- printf("Usage %s [extent] [beta] [num_iter] {[start]}\n", argv[0]);
- return -1;
- }
- // From the command line: set extent, beta and num_iter
- extent = atoi(argv[1]);
- beta = atof(argv[2]);
- num_iter= atoi(argv[3]);
- // Set the starting condition (to 0 if not set on command line)
- int start_cond = 0;
- if (argc==5)
- start_cond=atoi(argv[4]);
- // Check if all parameters are reasonble:
- if ((beta<0.0)||(beta>3.0)||(num_iter<1)) {
- printf("beta not in [0,3] (is: %f) or iterations to low (is %d)\n",
- beta, num_iter); return -1; }
- if ((extent<2)||(extent>128)||(extent%2)) {
- printf("The extent in not in [2,128] or not an even number (is %d)\n",
- extent); return -1; }
- if ((start_cond<0)||(start_cond>3)) {
- printf("Start condition: 0-random spins, 1-all up, 2-all down (is %d)\n",
- start_cond); return -1; }
- // Output file name gets the form "data/08/0.5160.dat" for extent=8,beta=0.516
- char out_file_name[64];
- mkdir("data", 0775);
- sprintf(out_file_name, "data/%02d", extent);
- mkdir(out_file_name, 0775);
- sprintf(out_file_name, "data/%02d/%.5f.dat", extent, beta);
- // Open the output file to write the results to
- FILE* outfile;
- outfile = fopen(out_file_name, "w");
- if (outfile == NULL) {
- printf("Output file could not be opened. Aborting.\n");
- return -1;
- }
- // Write some information to the output file
- fprintf(outfile, "## Results of the 2d Ising simulation\n");
- fprintf(outfile, "# beta= %1.4f\n# size= %2d\n# iter= %6d\n#start= %1d\n",
- beta, extent, num_iter, start_cond);
- // Initialize the spin array
- initialize(start_cond);
- printf("# Initialized, beta=%f, size=%2d, #iterations=%5d\n",
- beta, extent, num_iter);
- // The main update and measurement loop
- int i,j;
- for (i=0; i<num_iter; i++) {
- // Do 10 sweeps
- for (j=0;j<10;j++)
- sweep();
- // Run a measurement
- measure(outfile, i);
- // Output some status
- if ((i%1000)==0) {
- printf(".");
- fflush(stdout);
- }
- }
- printf(" done. \n");
- fclose(outfile);
- return 0;
- }
- #endif
- // Implementation of helper functions
- #ifndef _TEST_RUN
- double rnd() { return (drand48()); }
- int rnd_spin() { return ((rnd()<0.5)?(-1):(+1)); }
- void check_nn(int x, int y, int cur, int lns, int rns, int uns, int dns) { return; }
- #endif
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement