Trisurf Monte Carlo simulator
Samo Penic
2013-11-24 439a4e0f10dd120b563039012d07eb5f1b3c8ded
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];
50f10f 71     ts_uint i,j; //b;
SP 72 //    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 73     for(i=thID*partition_size;i<end;i++){
d7639a 74         rnvec[0]=drand48();
SP 75         rnvec[1]=drand48();
76         rnvec[2]=drand48();
c78004 77 //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 78
SP 79     /**** THREAD IS POTENTIALLY LOCKED ******/
439a4e 80     pthread_mutex_lock(&vesicle->mutex->vtx_taint); //taint if no other process is tainting or wait until you can taint
4e56fe 81     //    ts_fprintf(stdout, "thID=%d:: Tainting vertex %d, level=%d. Waiting....\n",thID, i, vesicle->vlist->vtx[i]->locked);
50f10f 82
c9884c 83     while(vertex_tainted(vesicle->vlist->vtx[i],1,1)){
c78004 84         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 85         for(j=0; j<vesicle->vlist->vtx[i]->neigh_no; j++){
c78004 86             ts_fprintf(stdout, "%d(a=%d) ", vesicle->vlist->vtx[i]->neigh[j]->idx,  vesicle->vlist->vtx[i]->neigh[j]->locked);
50f10f 87         }
c78004 88         ts_fprintf(stdout, ".\n"); 
200815 89     }
c78004 90     
439a4e 91     pthread_mutex_lock(&vesicle->mutex->vtx_untaint);
c78004 92     vertex_taint(vesicle->vlist->vtx[i],1);
439a4e 93     pthread_mutex_unlock(&vesicle->mutex->vtx_untaint);
50f10f 94
439a4e 95     pthread_mutex_unlock(&vesicle->mutex->vtx_taint);
4e56fe 96     /**** THREAD IS RELEASING MUTEX RESOURCES ******/
SP 97
c9884c 98         single_verticle_timestep(vesicle,vesicle->vlist->vtx[i],rnvec);
439a4e 99     pthread_mutex_lock(&vesicle->mutex->vtx_untaint);
c78004 100     vertex_untaint(vesicle->vlist->vtx[i],1);
439a4e 101     pthread_mutex_unlock(&vesicle->mutex->vtx_untaint);
200815 102 //        ts_fprintf(stdout, "Vertex %d should be untainted, level=%d.\n", i, vesicle->vlist->vtx[i]->locked);
d7639a 103     }
SP 104
b7c9b3 105 //    ts_int cnt=0;
c9884c 106 /*
142a67 107     for(i=0;i<vesicle->vlist->n;i++){
SP 108     b=rand() % vesicle->blist->n;
d7639a 109         //find a bond and return a pointer to a bond...
SP 110         //call single_bondflip_timestep...
142a67 111         retval=single_bondflip_timestep(vesicle,vesicle->blist->bond[b],rnvec);
b7c9b3 112 //    if(retval==TS_SUCCESS) cnt++;        
d7639a 113     } 
b7c9b3 114 //    printf("Bondflip success rate in one sweep: %d/%d=%e\n", cnt,vesicle->blist->n,(double)cnt/(double)vesicle->blist->n);
c9884c 115 */
SP 116 /*    if(retval);
d7639a 117     return TS_SUCCESS;
c9884c 118 */
SP 119      pthread_exit(0); /* exit */
d7639a 120 }
SP 121
122
123