Trisurf Monte Carlo simulator
view | history | commit | commitdiff
Samo Penic
2013-11-27 0652ee2715994b48a9bbfd1132bf0b65ad206289
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);
439a4e 25     pthread_mutex_init(&(vesicle->mutex->vtx_taint),NULL);
SP 26     pthread_mutex_init(&(vesicle->mutex->vtx_untaint),NULL);
27     pthread_mutex_init(&(vesicle->mutex->cell_modify),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     }
4e56fe 34
SP 35      for(i=0;i<inititer+iterations;i++){
d7a113 36         for(j=0;j<mcsweeps;j++){
50f10f 37             for(k=0;k<vesicle->threads;k++){
SP 38                 pthread_create(&tid[k], NULL, single_timestep, (void *)&data[k]);
39             }
40             for(k=0;k<vesicle->threads;k++){
41                     pthread_join(tid[k],NULL);
42             }
43         //    single_timestep(vesicle);
d7a113 44         }
SP 45         centermass(vesicle);
46         cell_occupation(vesicle);
47         if(i>inititer){
48             write_vertex_xml_file(vesicle,i-inititer);
49         }
50     }
4e56fe 51
439a4e 52     pthread_mutex_destroy(&(vesicle->mutex->vtx_taint));
SP 53     pthread_mutex_destroy(&(vesicle->mutex->vtx_untaint));
4e56fe 54     pthread_exit(NULL);
d7a113 55     return TS_SUCCESS;
SP 56 }
d7639a 57
c9884c 58 void * single_timestep(void *thread_data){
SP 59     thdata *data;
60     data=(thdata *)thread_data;
61     ts_vesicle *vesicle=data->vesicle;
62     ts_uint thID=data->thread_no;
50f10f 63     ts_uint partition_size=(ts_uint)(vesicle->vlist->n/data->threads);
c9884c 64     ts_uint end;
50f10f 65     if(thID==data->threads-1){
c9884c 66         end=vesicle->vlist->n;
SP 67     } else {
68         end=(thID+1)*partition_size;
69     }
d7639a 70     ts_double rnvec[3];
0652ee 71     ts_uint i,j,b;
SP 72     ts_vertex *vtx;
73     //ts_bool retval;    
50f10f 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 ******/
439a4e 82     pthread_mutex_lock(&vesicle->mutex->vtx_taint); //taint if no other process is tainting or wait until you can taint
4e56fe 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     
439a4e 93     pthread_mutex_lock(&vesicle->mutex->vtx_untaint);
c78004 94     vertex_taint(vesicle->vlist->vtx[i],1);
439a4e 95     pthread_mutex_unlock(&vesicle->mutex->vtx_untaint);
50f10f 96
439a4e 97     pthread_mutex_unlock(&vesicle->mutex->vtx_taint);
4e56fe 98     /**** THREAD IS RELEASING MUTEX RESOURCES ******/
SP 99
c9884c 100         single_verticle_timestep(vesicle,vesicle->vlist->vtx[i],rnvec);
439a4e 101     pthread_mutex_lock(&vesicle->mutex->vtx_untaint);
c78004 102     vertex_untaint(vesicle->vlist->vtx[i],1);
439a4e 103     pthread_mutex_unlock(&vesicle->mutex->vtx_untaint);
d7639a 104     }
SP 105
b7c9b3 106 //    ts_int cnt=0;
0652ee 107     
SP 108     for(i=0;i<(ts_uint)((ts_double)(vesicle->vlist->n)/(ts_double)vesicle->threads);i++){
142a67 109     b=rand() % vesicle->blist->n;
d7639a 110         //find a bond and return a pointer to a bond...
SP 111         //call single_bondflip_timestep...
0652ee 112     vtx=vesicle->blist->bond[b]->vtx1;
SP 113
114     /**** THREAD IS POTENTIALLY LOCKED ******/
115     pthread_mutex_lock(&vesicle->mutex->vtx_taint);
116             while(vertex_tainted(vtx,2,1)){
117         ts_fprintf(stdout, "BF:: thID=%d Vertex %d is tainted. Amount=%d. BU(H)TELJ!\n       neigh. vtxs: ",thID, i, *(vesicle->vlist->vtx[i]->locked));
118         for(j=0; j<vesicle->vlist->vtx[i]->neigh_no; j++){
119             ts_fprintf(stdout, "%d(a=%d) ", vesicle->vlist->vtx[i]->neigh[j]->idx,  *(vesicle->vlist->vtx[i]->neigh[j]->locked));
120         }
121         ts_fprintf(stdout, ".\n"); 
122     }
123     
124     pthread_mutex_lock(&vesicle->mutex->vtx_untaint);
125     vertex_taint(vtx,2);
126     pthread_mutex_unlock(&vesicle->mutex->vtx_untaint);
127
128     pthread_mutex_unlock(&vesicle->mutex->vtx_taint);
129     /**** THREAD IS RELEASING MUTEX RESOURCES ******/
130
131     //retval=
132     single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
133     
134     pthread_mutex_lock(&vesicle->mutex->vtx_untaint);
135     vertex_untaint(vtx,2);
136     pthread_mutex_unlock(&vesicle->mutex->vtx_untaint);
137
138
b7c9b3 139 //    if(retval==TS_SUCCESS) cnt++;        
d7639a 140     } 
b7c9b3 141 //    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n);
0652ee 142
c9884c 143 /*    if(retval);
d7639a 144     return TS_SUCCESS;
c9884c 145 */
SP 146      pthread_exit(0); /* exit */
d7639a 147 }
SP 148
149
150