pthredas working now (by partitioning atoms into count/#THREADS blocks)
[physik/posic.git] / moldyn.c
index 6da2b95..ed30815 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
 #include <omp.h>
 #endif
 
+#if defined PTHREADS || defined VISUAL_THREAD
+#include <pthread.h>
+#endif
+
 #include "moldyn.h"
 #include "report/report.h"
 
 #undef PSE_NAME
 #undef PSE_COL
 
+#ifdef PTHREADS
+/* global mutexes */
+pthread_mutex_t *amutex;
+pthread_mutex_t emutex;
+#endif
+
 /*
  * the moldyn functions
  */
@@ -62,13 +72,27 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
        rand_init(&(moldyn->random),NULL,1);
        moldyn->random.status|=RAND_STAT_VERBOSE;
 
+#ifdef PTHREADS
+       pthread_mutex_init(&emutex,NULL);
+#endif
+
        return 0;
 }
 
 int moldyn_shutdown(t_moldyn *moldyn) {
 
+#ifdef PTHREADS
+       int i;
+#endif
+
        printf("[moldyn] shutdown\n");
 
+#ifdef PTHREADS
+       for(i=0;i<moldyn->count;i++)
+               pthread_mutex_destroy(&(amutex[i]));
+       pthread_mutex_destroy(&emutex);
+#endif
+
        moldyn_log_shutdown(moldyn);
        link_cell_shutdown(moldyn);
        rand_close(&(moldyn->random));
@@ -496,6 +520,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        void *ptr;
        t_atom *atom;
        char name[16];
+#ifdef PTHREADS
+       pthread_mutex_t *mutex;
+#endif
 
        new=a*b*c;
        count=moldyn->count;
@@ -518,6 +545,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        moldyn->atom=ptr;
        atom=&(moldyn->atom[count]);
 
+#ifdef PTHREADS
+       ptr=realloc(amutex,(count+new)*sizeof(pthread_mutex_t));
+       if(!ptr) {
+               perror("[moldyn] mutex realloc (add atom)");
+               return -1;
+       }
+       amutex=ptr;
+       mutex=&(amutex[count]);
+#endif
+
        /* no atoms on the boundaries (only reason: it looks better!) */
        if(!origin) {
                orig.x=0.5*lc;
@@ -575,6 +612,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                atom[ret].tag=count+ret;
                check_per_bound(moldyn,&(atom[ret].r));
                atom[ret].r_0=atom[ret].r;
+#ifdef PTHREADS
+               pthread_mutex_init(&(mutex[ret]),NULL);
+#endif
        }
 
        /* update total system mass */
@@ -609,6 +649,16 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        moldyn->lc.subcell->list=ptr;
 #endif
 
+#ifdef PTHREADS
+       ptr=realloc(amutex,(count+1)*sizeof(pthread_mutex_t));
+       if(!ptr) {
+               perror("[moldyn] mutex realloc (add atom)");
+               return -1;
+       }
+       amutex=ptr;
+       pthread_mutex_init(&(amutex[count]),NULL);
+#endif
+
        atom=moldyn->atom;
 
        /* initialize new atom */
@@ -1620,6 +1670,17 @@ int moldyn_integrate(t_moldyn *moldyn) {
        struct timeval t1,t2;
        //double tp;
 
+#ifdef VISUAL_THREAD
+       u8 first,change;
+       pthread_t io_thread;
+       int ret;
+       t_moldyn md_copy;
+       t_atom *atom_copy;
+
+       first=1;
+       change=0;
+#endif
+
        sched=&(moldyn->schedule);
        atom=moldyn->atom;
 
@@ -1668,6 +1729,16 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* debugging, ignore */
        moldyn->debug=0;
 
+       /* zero & init moldyn copy */
+#ifdef VISUAL_THREAD
+       memset(&md_copy,0,sizeof(t_moldyn));
+       atom_copy=malloc(moldyn->count*sizeof(t_atom));
+       if(atom_copy==NULL) {
+               perror("[moldyn] malloc atom copy (init)");
+               return -1;
+       }
+#endif
+
        /* tell the world */
        printf("[moldyn] integration start, go get a coffee ...\n");
 
@@ -1752,7 +1823,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
                if(s) {
                        if(!(moldyn->total_steps%s)) {
-                               snprintf(dir,128,"%s/s-%07.f.save",
+                               snprintf(dir,128,"%s/s-%08.f.save",
                                         moldyn->vlsdir,moldyn->time);
                                fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
                                        S_IRUSR|S_IWUSR);
@@ -1767,16 +1838,45 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
                if(a) {
                        if(!(moldyn->total_steps%a)) {
+#ifdef VISUAL_THREAD
+       /* check whether thread has not terminated yet */
+       if(!first) {
+               ret=pthread_join(io_thread,NULL);
+       }
+       first=0;
+       /* prepare and start thread */
+       if(moldyn->count!=md_copy.count) {
+               free(atom_copy);
+               change=1;
+       }
+       memcpy(&md_copy,moldyn,sizeof(t_moldyn));
+       if(change) {
+               atom_copy=malloc(moldyn->count*sizeof(t_atom));
+               if(atom_copy==NULL) {
+                       perror("[moldyn] malloc atom copy (change)");
+                       return -1;
+               }
+       }
+       md_copy.atom=atom_copy;
+       memcpy(atom_copy,moldyn->atom,moldyn->count*sizeof(t_atom));
+       change=0;
+       ret=pthread_create(&io_thread,NULL,visual_atoms,&md_copy);
+       if(ret) {
+               perror("[moldyn] create visual atoms thread\n");
+               return -1;
+       }
+#else
                                visual_atoms(moldyn);
+#endif
                        }
                }
 
                /* display progress */
-               //if(!(moldyn->total_steps%10)) {
+               if(!(i%10)) {
                        /* get current time */
                        gettimeofday(&t2,NULL);
 
-printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
        moldyn->p/BAR,moldyn->p_avg/BAR,
@@ -1788,7 +1888,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
 
                        /* copy over time */
                        t1=t2;
-               //}
+               }
 
                /* increase absolute time */
                moldyn->time+=moldyn->tau;
@@ -2117,6 +2217,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                                  jtom,
                                                                  ktom,
                                                                  bc_ik|bc_ij);
+
 #ifdef STATIC_LISTS
                                        }
 #elif LOWMEM_LISTS
@@ -2515,6 +2616,84 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 
 }
 
+/*
+ * function to find neighboured atoms
+ */
+
+int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom,
+                       int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom,
+                                     void *data,u8 bc)) {
+
+       t_linkcell *lc;
+#ifdef STATIC_LISTS
+       int *neighbour[27];
+       int p;
+#elif LOWMEM_LISTS
+       int neighbour[27];
+       int p;
+#else
+       t_list neighbour[27];
+       t_list *this;
+#endif
+       u8 bc;
+       t_atom *natom;
+       int j;
+
+       lc=&(moldyn->lc);
+       
+       /* neighbour indexing */
+       link_cell_neighbour_index(moldyn,
+                                 (atom->r.x+moldyn->dim.x/2)/lc->x,
+                                 (atom->r.y+moldyn->dim.y/2)/lc->x,
+                                 (atom->r.z+moldyn->dim.z/2)/lc->x,
+                                 neighbour);
+
+       for(j=0;j<27;j++) {
+
+               bc=(j<lc->dnlc)?0:1;
+
+#ifdef STATIC_LISTS
+               p=0;
+
+               while(neighbour[j][p]!=-1) {
+
+                       natom=&(moldyn->atom[neighbour[j][p]]);
+                       p++;
+#elif LOWMEM_LISTS
+               p=neighbour[j];
+
+               while(p!=-1) {
+
+                       natom=&(moldyn->atom[p]);
+                       p=lc->subcell->list[p];
+#else
+               this=&(neighbour[j]);
+               list_reset_f(this);
+
+               if(this->start==NULL)
+                       continue;
+
+               do {
+
+                       natom=this->current->data;
+#endif
+
+                       /* process bond */
+                       process(moldyn,atom,natom,data,bc);
+
+#ifdef STATIC_LISTS
+               }
+#elif LOWMEM_LISTS
+               }
+#else
+               } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
+       }
+
+       return 0;
+
+}
+
 /*
  * post processing functions
  */
@@ -2841,7 +3020,11 @@ int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
        return 0;
 }
 
+#ifdef VISUAL_THREAD
+void *visual_atoms(void *ptr) {
+#else
 int visual_atoms(t_moldyn *moldyn) {
+#endif
 
        int i;
        char file[128+64];
@@ -2850,6 +3033,11 @@ int visual_atoms(t_moldyn *moldyn) {
        t_visual *v;
        t_atom *atom;
        t_vb vb;
+#ifdef VISUAL_THREAD
+       t_moldyn *moldyn;
+
+       moldyn=ptr;
+#endif
 
        v=&(moldyn->vis);
        dim.x=v->dim.x;
@@ -2859,17 +3047,19 @@ int visual_atoms(t_moldyn *moldyn) {
 
        help=(dim.x+dim.y);
 
-       sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
+       sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time);
        vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
        if(vb.fd<0) {
                perror("open visual save file fd");
+#ifndef VISUAL_THREAD
                return -1;
+#endif
        }
 
        /* write the actual data file */
 
        // povray header
-       dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
+       dprintf(vb.fd,"# [P] %d %08.f <%f,%f,%f>\n",
                moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
 
        // atomic configuration
@@ -2883,7 +3073,9 @@ int visual_atoms(t_moldyn *moldyn) {
                                                    atom[i].ekin);
        
        // bonds between atoms
+#ifndef VISUAL_THREAD
        process_2b_bonds(moldyn,&vb,visual_bonds_process);
+#endif
        
        // boundaries
        if(dim.x) {
@@ -2929,8 +3121,15 @@ int visual_atoms(t_moldyn *moldyn) {
 
        close(vb.fd);
 
+#ifdef VISUAL_THREAD
+       pthread_exit(NULL);
+
+}
+#else
+
        return 0;
 }
+#endif
 
 /*
  * fpu cntrol functions