Trisurf Monte Carlo simulator
Samo Penic
2013-11-23 c78004bcfacb0dc6399647027b3b1acc013db153
commit | author | age
d7639a 1 #include<stdlib.h>
SP 2 #include<stdio.h>
aec47d 3 #include<math.h>
c9884c 4 #include<pthread.h>
c78004 5 //#include<unistd.h> //usleep requires it
aec47d 6 //#include "io.h"
SP 7 #include "general.h"
8 #include "timestep.h"
9 #include "vertexmove.h"
30ee9c 10 #include "bondflip.h"
d7a113 11 #include "frame.h"
200815 12 #include "vertex.h"
d7a113 13 #include "io.h"
c9884c 14
d7a113 15 ts_bool run_simulation(ts_vesicle *vesicle, ts_uint mcsweeps, ts_uint inititer, ts_uint iterations){
50f10f 16     ts_uint i, j,k ;
d7a113 17
SP 18     centermass(vesicle);
19     cell_occupation(vesicle);
20     ts_fprintf(stdout, "Starting simulation (first %d x %d MC sweeps will not be recorded on disk)\n", inititer, mcsweeps);
50f10f 21     ts_fprintf(stdout, "Starting %d threads\n",vesicle->threads);
4e56fe 22
50f10f 23     pthread_t *tid=(pthread_t *)malloc(sizeof(pthread_t)*vesicle->threads);
SP 24     thdata *data=(thdata *)malloc(sizeof(thdata)*vesicle->threads);
25     pthread_mutex_t mutex_taint, mutex_untaint;
26     pthread_mutex_init(&mutex_taint,NULL);
27     pthread_mutex_init(&mutex_untaint,NULL);
4e56fe 28
50f10f 29     for(i=0;i<vesicle->threads;i++){
SP 30         data[i].thread_no=i;
31         data[i].vesicle=vesicle;
32         data[i].threads=vesicle->threads;
33         data[i].vtx_tainting_lock=&mutex_taint;
34         data[i].vtx_untainting_lock=&mutex_untaint;
35     }
4e56fe 36
SP 37      for(i=0;i<inititer+iterations;i++){
d7a113 38         for(j=0;j<mcsweeps;j++){
50f10f 39             for(k=0;k<vesicle->threads;k++){
SP 40                 pthread_create(&tid[k], NULL, single_timestep, (void *)&data[k]);
41             }
42             for(k=0;k<vesicle->threads;k++){
43                     pthread_join(tid[k],NULL);
44             }
45         //    single_timestep(vesicle);
d7a113 46         }
SP 47         centermass(vesicle);
48         cell_occupation(vesicle);
49         if(i>inititer){
50             write_vertex_xml_file(vesicle,i-inititer);
51         }
52     }
4e56fe 53
50f10f 54     pthread_mutex_destroy(&mutex_taint);
SP 55     pthread_mutex_destroy(&mutex_untaint);
4e56fe 56     pthread_exit(NULL);
d7a113 57     return TS_SUCCESS;
SP 58 }
d7639a 59
c9884c 60 void * single_timestep(void *thread_data){
SP 61     thdata *data;
62     data=(thdata *)thread_data;
63     ts_vesicle *vesicle=data->vesicle;
64     ts_uint thID=data->thread_no;
50f10f 65     ts_uint partition_size=(ts_uint)(vesicle->vlist->n/data->threads);
c9884c 66     ts_uint end;
50f10f 67     if(thID==data->threads-1){
c9884c 68         end=vesicle->vlist->n;
SP 69     } else {
70         end=(thID+1)*partition_size;
71     }
d7639a 72     ts_double rnvec[3];
50f10f 73     ts_uint i,j; //b;
SP 74 //    ts_fprintf(stdout,"Thread thID=%d report for duty. My vtxes are in range from %d to %d.\n",thID,thID*partition_size, end-1);
c9884c 75     for(i=thID*partition_size;i<end;i++){
d7639a 76         rnvec[0]=drand48();
SP 77         rnvec[1]=drand48();
78         rnvec[2]=drand48();
c78004 79 //TODO: you only need to taint (lock) vertices, that could be shared with another partition. You can leave out the rest. How can you test for that? Well, if any of vtx neighbours are in other partition, then you need to lock vertex and that neighbor. Otherwise not!
4e56fe 80
SP 81     /**** THREAD IS POTENTIALLY LOCKED ******/
82     pthread_mutex_lock(data->vtx_tainting_lock); //taint if no other process is tainting or wait until you can taint
83     //    ts_fprintf(stdout, "thID=%d:: Tainting vertex %d, level=%d. Waiting....\n",thID, i, vesicle->vlist->vtx[i]->locked);
50f10f 84
c9884c 85     while(vertex_tainted(vesicle->vlist->vtx[i],1,1)){
c78004 86         ts_fprintf(stdout, "thID=%d:: Vertex %d is tainted, so I cannot try vertexmove. Amount=%d. BU(H)TELJ!\n       neigh. vtxs: ",thID, i, vesicle->vlist->vtx[i]->locked);
50f10f 87         for(j=0; j<vesicle->vlist->vtx[i]->neigh_no; j++){
c78004 88             ts_fprintf(stdout, "%d(a=%d) ", vesicle->vlist->vtx[i]->neigh[j]->idx,  vesicle->vlist->vtx[i]->neigh[j]->locked);
50f10f 89         }
c78004 90         ts_fprintf(stdout, ".\n"); 
200815 91     }
c78004 92     
50f10f 93     pthread_mutex_lock(data->vtx_untainting_lock);
c78004 94     vertex_taint(vesicle->vlist->vtx[i],1);
50f10f 95     pthread_mutex_unlock(data->vtx_untainting_lock);
SP 96
4e56fe 97     pthread_mutex_unlock(data->vtx_tainting_lock);
SP 98     /**** THREAD IS RELEASING MUTEX RESOURCES ******/
99
c9884c 100         single_verticle_timestep(vesicle,vesicle->vlist->vtx[i],rnvec);
50f10f 101     pthread_mutex_lock(data->vtx_untainting_lock);
c78004 102     vertex_untaint(vesicle->vlist->vtx[i],1);
50f10f 103     pthread_mutex_unlock(data->vtx_untainting_lock);
200815 104 //        ts_fprintf(stdout, "Vertex %d should be untainted, level=%d.\n", i, vesicle->vlist->vtx[i]->locked);
d7639a 105     }
SP 106
b7c9b3 107 //    ts_int cnt=0;
c9884c 108 /*
142a67 109     for(i=0;i<vesicle->vlist->n;i++){
SP 110     b=rand() % vesicle->blist->n;
d7639a 111         //find a bond and return a pointer to a bond...
SP 112         //call single_bondflip_timestep...
142a67 113         retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
b7c9b3 114 //    if(retval==TS_SUCCESS) cnt++;        
d7639a 115     } 
b7c9b3 116 //    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n);
c9884c 117 */
SP 118 /*    if(retval);
d7639a 119     return TS_SUCCESS;
c9884c 120 */
SP 121      pthread_exit(0); /* exit */
d7639a 122 }
SP 123
124
125