Advertisement
Guest User

Untitled

a guest
Jul 26th, 2017
95
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 11.54 KB | None | 0 0
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <sys/stat.h>
  5.  
  6. // =======================================================================
  7. // Global variables and helper functions you can use
  8. // -----------------------------------------------------------------------
  9. signed char spins[128][128];    //!< The spin array: either -1 (down) or 1 (up)
  10. int   extent;           //!< The extent (size) of the lattice
  11. float beta;             //!< \beta
  12. int   num_iter;         //!< number of iterations
  13.  
  14. double rnd();       //!< returns a random number in [0,1]
  15. int rnd_spin();     //!< returns a random spin, either -1 or 1
  16.  
  17. //! Used by the test routine, checks if you set the next neighbours correctly
  18. void check_nn(int, int, int, int, int, int, int) ;
  19.  
  20.  
  21.  
  22. // =======================================================================
  23. // Declarations of all functions you need to edit
  24. // -----------------------------------------------------------------------
  25. void initialize(int);   //!< initialize spins
  26. void sweep();       //!< run an update sweep
  27. void measure();     //!< measure magnetisation
  28. int  calculate_spin(int, int, int, int, int); //!< return a new spin
  29.  
  30.  
  31.  
  32. // =======================================================================
  33. // A routine to initalize the spin array
  34. // -----------------------------------------------------------------------
  35. void initialize(int start_cond) {
  36.  
  37.     // Loop through the spin array
  38.     int x,y;
  39.     for(x=0; x<extent; x++)
  40.         for(y=0; y<extent; y++) {
  41.  
  42.             // /////////////////////////////////////////////////////////////////////////
  43.  
  44.             /* Task 1:  To warm up, something easy first: Depending on the starting
  45.              * condition (variable "start_cond"), the spins should be
  46.              * initialized:
  47.              *
  48.              * 0: random spins, 1: all spins up (1), 2: all spins down (-1)
  49.              *
  50.              * A loop over the lattice is alreaddy there, you only need
  51.              * to check the starting condition (e.g. with an if-statement)
  52.              * and set spin[x][y] accordingly. Use rnd_spin() to get random
  53.              * spins....
  54.              */
  55.  
  56.  
  57.             // your code here ...
  58.             if (start_cond == 0) {
  59.                 spins[x][y] = rnd_spin();
  60.             }
  61.             else if (start_cond == 1) {
  62.                 spins[x][y] = 1;
  63.             }
  64.             else if (start_cond == 2) {
  65.                 spins[x][y] = -1;
  66.             }
  67.  
  68.             // /////////////////////////////////////////////////////////////////////////
  69.  
  70.         }
  71. }
  72.  
  73.  
  74.  
  75. // =======================================================================
  76. // A routine to update the spin array
  77. // -----------------------------------------------------------------------
  78. void sweep() {
  79.  
  80.     // variables
  81.     int cur;            // the current spin
  82.     int lns, rns, uns, dns; // left, right, up and down neighbour to the spin
  83.  
  84.     int x,y,eo;         // loop variables
  85.     for (eo=0;eo<=1;eo++)       // even/odd splitting
  86.         for (x=0;x<extent;x++)      // loop x
  87.             for (y=0;y<extent;y++) {    // loop y
  88.  
  89.                 // Do nothing if the coordinate (x,y) we are currently
  90.                 // on does not match the even/odd variable
  91.                 if (((x+y)%2)!=eo) continue;
  92.  
  93.                 // get the current spin
  94.                 cur  = spins[x][y];
  95.  
  96.                 // /////////////////////////////////////////////////////////////////////////
  97.  
  98.                 /* Task 2:  This is inside a loop over the whole lattice. The variable "cur"
  99.                  * now contains the spin at x,y. You must fill the variables
  100.                  * "lns", "rns", "uns", "dns" with the correct neighbours:
  101.  
  102.                  * left: x-1, right: x+1, up: y-1, down: y+1
  103.                  *
  104.                  * lns = ....
  105.                  * rns = ....
  106.                  *
  107.                  * Remember that there are boundary conditions, so take care
  108.                  * for x,y=0 and x,y=(extent-1)
  109.                  */
  110.  
  111.  
  112.                 // your code here
  113.                 if (x > 0 && x < (extent-1) && y > 0 && y < (extent-1)) {
  114.                     lns = spins[x-1][y];
  115.                     rns = spins[x+1][y];
  116.                     uns = spins[x][y-1];
  117.                     dns = spins[x][y+1];
  118.                 }
  119.                 else if (y == 0 && x > 0 && x < (extent-1)) {
  120.                     lns = spins[x-1][y];
  121.                     rns = spins[x+1][y];
  122.                     dns = spins[x][y+1];
  123.                 }
  124.                 else if (y == (extent-1) && x > 0 && x < (extent-1)) {
  125.                     lns = spins[x-1][y];
  126.                     rns = spins[x+1][y];
  127.                     uns = spins[x][y-1];
  128.                 }
  129.                 else if (x == 0 && y > 0 && y < (extent-1)) {
  130.                     rns = spins[x+1][y];
  131.                     uns = spins[x][y-1];
  132.                     dns = spins[x][y+1];
  133.                 }
  134.                 else if (x == (extent-1) && y > 0 && y < (extent-1)) {
  135.                     lns = spins[x-1][y];
  136.                     uns = spins[x][y-1];
  137.                     dns = spins[x][y+1];
  138.                 }
  139.                 else if (x == 0 && y == 0) {
  140.                     rns = spins[x+1][y];
  141.                     dns = spins[x][y+1];
  142.                 }
  143.                 else if (x == 0 && y == (extent-1)) {
  144.                     rns = spins[x+1][y];
  145.                     uns = spins[x][y-1];
  146.                 }
  147.                 else if (x == (extent-1) && y == 0) {
  148.                     lns = spins[x-1][y];
  149.                     dns = spins[x][y+1];
  150.                 }
  151.                 else if (x == (extent-1) && y == (extent-1)) {
  152.                     lns = spins[x-1][y];
  153.                     uns = spins[x][y-1];
  154.                 }
  155.  
  156.                 // /////////////////////////////////////////////////////////////////////////
  157.  
  158.                 // if you compile the test, this checks if you set the neighbours right
  159.                 // if you compile the normal programm, this does nothing
  160.                 check_nn(x,y,eo,lns,rns,uns,dns);
  161.  
  162.                 // update the current spin (you will have to complete Task 3):
  163.                 spins[x][y] = calculate_spin(cur, lns, rns, uns, dns);
  164.             }
  165. }
  166.  
  167.  
  168. // =======================================================================
  169. // Calculate a new spin from its neighbours and beta
  170. // -----------------------------------------------------------------------
  171. int calculate_spin(int cur, int lns, int rns, int uns, int dns) {
  172.  
  173.     int new_spin;
  174.  
  175.     // ////////////////////////////////////////////////////////////////////////
  176.  
  177.     /* Task 3 : Now some physics:
  178.  
  179.      * Here you have to calculate a new spin. "cur" holds the spin, "lns",
  180.      * "rns", "uns" and "dns" its neighbours (o.k., you probably know that
  181.      * because you just wrote that part ....)
  182.      *
  183.      * As you (should) know, in the Ising model the spin flips if the flipped
  184.      * position is energetically favored or equal. If it energetically
  185.      * disfavored, there is still a (temperature-dependent) chance
  186.      * the spin will flip:
  187.      *
  188.      * P_flip = exp(-beta * (energy difference))
  189.      *
  190.      * You can (and should) use "rnd()" to generate a random number
  191.      * between 0 and 1 to implement this.
  192.      *
  193.      */
  194.  
  195.  
  196.     // your code here ...
  197.  
  198.  
  199.     // /////////////////////////////////////////////////////////////////////////
  200.  
  201.     return new_spin;
  202. }
  203.  
  204.  
  205. // =======================================================================
  206. // Measures and outputs energy and magnetisation to file
  207. // (it's already complete, you do not need to edit it)
  208. // -----------------------------------------------------------------------
  209. void measure(FILE* output, int iteration) {
  210.  
  211.     int x,y;
  212.  
  213.     double ene=0; // Store enery and magnetisation
  214.     double mag=0;
  215.  
  216.     // Energy via the Hamiltonian - no boundary-crossing terms,
  217.     // but on average they don't matter
  218.     for(y=0;y<extent-1;y++)
  219.         for(x=0;x<extent-1;x++)
  220.             ene += 2*spins[x][y] * ( spins[x][y+1] + spins[x+1][y] );
  221.  
  222.     // Magnetisation as sum of spin directions
  223.     for(y=0;y<extent;y++)
  224.         for(x=0;x<extent;x++)
  225.             mag += spins[x][y];
  226.  
  227.     // Both have to be divided, so they are "per site" values
  228.     ene /= (extent-1)*(extent-1);
  229.     mag /= extent*extent;
  230.  
  231.     // Write the result to the output file
  232.     fprintf(output, "%6d\t%e\t%e\n", iteration, ene, fabs(mag));
  233. }
  234.  
  235.  
  236.  
  237. /* ****************** END ************************************************
  238.  *  
  239.  * Below here, the main routine that controls the programm and the
  240.  * helper functions are implemented. You do not need to change anything
  241.  * here, you only need to read these if you are interested.
  242.  *
  243.  * ********************************************************************** */
  244.  
  245. #ifndef _TEST_RUN //Used for test-ising, don't modify these "#ifndef" things
  246.  
  247. //! The main routine controlling program flow
  248. int main(int argc, char** argv) {
  249.  
  250.     // Check the number of command line arguments
  251.     if ((argc!=4)&&(argc!=5)) {
  252.         printf("Usage %s [extent] [beta] [num_iter] {[start]}\n", argv[0]);
  253.         return -1;
  254.     }
  255.  
  256.     // From the command line: set extent, beta and num_iter
  257.     extent  = atoi(argv[1]);
  258.     beta    = atof(argv[2]);
  259.     num_iter= atoi(argv[3]);
  260.  
  261.     // Set the starting condition (to 0 if not set on command line)
  262.     int start_cond = 0;
  263.     if (argc==5)
  264.         start_cond=atoi(argv[4]);
  265.  
  266.     // Check if all parameters are reasonble:
  267.     if ((beta<0.0)||(beta>3.0)||(num_iter<1)) {
  268.         printf("beta not in [0,3] (is: %f) or iterations to low (is %d)\n",
  269.                 beta, num_iter); return -1; }
  270.  
  271.     if ((extent<2)||(extent>128)||(extent%2)) {
  272.         printf("The extent in not in [2,128] or not an even number (is %d)\n",
  273.                 extent); return -1; }
  274.  
  275.     if ((start_cond<0)||(start_cond>3)) {
  276.         printf("Start condition: 0-random spins, 1-all up, 2-all down (is %d)\n",
  277.                 start_cond); return -1; }
  278.  
  279.     // Output file name gets the form "data/08/0.5160.dat" for extent=8,beta=0.516
  280.     char out_file_name[64];
  281.     mkdir("data", 0775);
  282.     sprintf(out_file_name, "data/%02d", extent);
  283.     mkdir(out_file_name, 0775);
  284.     sprintf(out_file_name, "data/%02d/%.5f.dat", extent, beta);
  285.  
  286.     // Open the output file to write the results to
  287.     FILE* outfile;
  288.     outfile = fopen(out_file_name, "w");
  289.     if (outfile == NULL) {
  290.         printf("Output file could not be opened. Aborting.\n");
  291.         return -1;
  292.     }
  293.  
  294.     // Write some information to the output file
  295.     fprintf(outfile, "## Results of the 2d Ising simulation\n");
  296.     fprintf(outfile, "# beta= %1.4f\n# size= %2d\n# iter= %6d\n#start= %1d\n",
  297.             beta, extent, num_iter, start_cond);
  298.  
  299.     // Initialize the spin array
  300.     initialize(start_cond);
  301.     printf("# Initialized, beta=%f, size=%2d, #iterations=%5d\n",
  302.             beta, extent, num_iter);
  303.  
  304.     // The main update and measurement loop
  305.     int i,j;
  306.     for (i=0; i<num_iter; i++) {
  307.         // Do 10 sweeps
  308.         for (j=0;j<10;j++)
  309.             sweep();
  310.         // Run a measurement
  311.         measure(outfile, i);
  312.         // Output some status
  313.         if ((i%1000)==0) {
  314.             printf(".");
  315.             fflush(stdout);
  316.         }
  317.  
  318.     }
  319.  
  320.     printf(" done. \n");    
  321.     fclose(outfile);
  322.     return 0;
  323. }
  324. #endif
  325.  
  326. // Implementation of helper functions
  327. #ifndef _TEST_RUN
  328. double rnd()    {  return (drand48());  }
  329. int rnd_spin()  {  return ((rnd()<0.5)?(-1):(+1));  }
  330. void check_nn(int x, int y, int cur, int lns, int rns, int uns, int dns) {  return; }
  331. #endif
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement