X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=8451c4436f43db2213bec83677c359fb4062f133;hb=787848dc9bf3daa2ccb76296de904688cc504a45;hp=e936b11d0ce5a977330d9989bc7c156cfc25e014;hpb=45f9190fc7489f3f9ff8f0d917945b1d5bc9d584;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index e936b11..8451c44 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -23,7 +23,7 @@
 #include <omp.h>
 #endif
 
-#ifdef PTHREADS
+#if defined PTHREADS || defined VISUAL_THREAD
 #include <pthread.h>
 #endif
 
@@ -47,6 +47,12 @@
 #undef PSE_NAME
 #undef PSE_COL
 
+#ifdef PTHREADS
+/* global mutexes */
+pthread_mutex_t *amutex;
+pthread_mutex_t emutex;
+#endif
+
 /*
  * the moldyn functions
  */
@@ -66,13 +72,28 @@ 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]));
+	free(amutex);
+	pthread_mutex_destroy(&emutex);
+#endif
+
 	moldyn_log_shutdown(moldyn);
 	link_cell_shutdown(moldyn);
 	rand_close(&(moldyn->random));
@@ -492,7 +513,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  */
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
+                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
+                   u8 p_type,t_part_vals *p_vals) {
 
 	int new,count;
 	int ret;
@@ -500,6 +522,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;
@@ -522,6 +547,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;
@@ -537,21 +572,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	switch(type) {
 		case CUBIC:
 			set_nn_dist(moldyn,lc);
-			ret=cubic_init(a,b,c,lc,atom,&orig);
+			ret=cubic_init(a,b,c,lc,atom,&orig,p_type,p_vals);
 			strcpy(name,"cubic");
 			break;
 		case FCC:
 			if(!origin)
 				v3_scale(&orig,&orig,0.5);
 			set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
-			ret=fcc_init(a,b,c,lc,atom,&orig);
+			ret=fcc_init(a,b,c,lc,atom,&orig,p_type,p_vals);
 			strcpy(name,"fcc");
 			break;
 		case DIAMOND:
 			if(!origin)
 				v3_scale(&orig,&orig,0.25);
 			set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
-			ret=diamond_init(a,b,c,lc,atom,&orig);
+			ret=diamond_init(a,b,c,lc,atom,&orig,p_type,p_vals);
 			strcpy(name,"diamond");
 			break;
 		default:
@@ -579,6 +614,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 */
@@ -613,6 +651,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 */
@@ -663,12 +711,14 @@ int del_atom(t_moldyn *moldyn,int tag) {
 }
 
 /* cubic init */
-int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+               u8 p_type,t_part_vals *p_vals) {
 
 	int count;
 	t_3dvec r;
 	int i,j,k;
 	t_3dvec o;
+	t_3dvec dist;
 
 	count=0;
 	if(origin)
@@ -676,14 +726,39 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	else
 		v3_zero(&o);
 
+	/* shift partition values */
+	if(p_type) {
+		p_vals->p.x+=(a*lc)/2.0;
+		p_vals->p.y+=(b*lc)/2.0;
+		p_vals->p.z+=(c*lc)/2.0;
+	}
+
 	r.x=o.x;
 	for(i=0;i<a;i++) {
 		r.y=o.y;
 		for(j=0;j<b;j++) {
 			r.z=o.z;
 			for(k=0;k<c;k++) {
+				switch(p_type) {
+					case PART_INSIDE_R:
+						v3_sub(&dist,&r,&(p_vals->p));
+			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
 				v3_copy(&(atom[count].r),&r);
 				count+=1;
+			}
+						break;
+					case PART_OUTSIDE_R:
+						v3_sub(&dist,&r,&(p_vals->p));
+			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+			}
+						break;
+					default:	
+						v3_copy(&(atom[count].r),&r);
+						count+=1;
+						break;
+				}
 				r.z+=lc;
 			}
 			r.y+=lc;
@@ -701,12 +776,14 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 }
 
 /* fcc lattice init */
-int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+             u8 p_type,t_part_vals *p_vals) {
 
 	int count;
 	int i,j,k,l;
 	t_3dvec o,r,n;
 	t_3dvec basis[3];
+	t_3dvec dist;
 
 	count=0;
 	if(origin)
@@ -723,6 +800,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	basis[2].y=0.5*lc;
 	basis[2].z=0.5*lc;
 
+	/* shift partition values */
+	if(p_type) {
+		p_vals->p.x+=(a*lc)/2.0;
+		p_vals->p.y+=(b*lc)/2.0;
+		p_vals->p.z+=(c*lc)/2.0;
+	}
+
 	/* fill up the room */
 	r.x=o.x;
 	for(i=0;i<a;i++) {
@@ -731,15 +815,51 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 			r.z=o.z;
 			for(k=0;k<c;k++) {
 				/* first atom */
+				switch(p_type) {
+					case PART_INSIDE_R:
+						v3_sub(&dist,&r,&(p_vals->p));
+			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
 				v3_copy(&(atom[count].r),&r);
 				count+=1;
-				r.z+=lc;
+			}
+						break;
+					case PART_OUTSIDE_R:
+						v3_sub(&dist,&r,&(p_vals->p));
+			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+			}
+						break;
+					default:
+						v3_copy(&(atom[count].r),&r);
+						count+=1;
+						break;
+				}
 				/* the three face centered atoms */
 				for(l=0;l<3;l++) {
 					v3_add(&n,&r,&basis[l]);
+					switch(p_type) {
+						case PART_INSIDE_R:
+			v3_sub(&dist,&n,&(p_vals->p));
+			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+			}
+							break;
+						case PART_OUTSIDE_R:
+			v3_sub(&dist,&n,&(p_vals->p));
+			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+			}
+							break;
+						default:
 					v3_copy(&(atom[count].r),&n);
 					count+=1;
+							break;
+					}
 				}
+				r.z+=lc;
 			}
 			r.y+=lc;
 		}
@@ -756,12 +876,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	return count;
 }
 
-int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+                 u8 p_type,t_part_vals *p_vals) {
 
 	int count;
 	t_3dvec o;
 
-	count=fcc_init(a,b,c,lc,atom,origin);
+	count=fcc_init(a,b,c,lc,atom,origin,p_type,p_vals);
 
 	o.x=0.25*lc;
 	o.y=0.25*lc;
@@ -769,7 +890,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
 	if(origin) v3_add(&o,&o,origin);
 
-	count+=fcc_init(a,b,c,lc,&atom[count],&o);
+	count+=fcc_init(a,b,c,lc,&atom[count],&o,p_type,p_vals);
 
 	return count;
 }
@@ -1624,7 +1745,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	struct timeval t1,t2;
 	//double tp;
 
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 	u8 first,change;
 	pthread_t io_thread;
 	int ret;
@@ -1684,7 +1805,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	moldyn->debug=0;
 
 	/* zero & init moldyn copy */
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 	memset(&md_copy,0,sizeof(t_moldyn));
 	atom_copy=malloc(moldyn->count*sizeof(t_atom));
 	if(atom_copy==NULL) {
@@ -1693,6 +1814,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	}
 #endif
 
+#ifdef PTHREADS
+	printf("##################\n");
+	printf("# USING PTHREADS #\n");
+	printf("##################\n");
+#endif
 	/* tell the world */
 	printf("[moldyn] integration start, go get a coffee ...\n");
 
@@ -1792,7 +1918,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 		if(a) {
 			if(!(moldyn->total_steps%a)) {
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 	/* check whether thread has not terminated yet */
 	if(!first) {
 		ret=pthread_join(io_thread,NULL);
@@ -1830,7 +1956,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			/* 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,
@@ -2171,6 +2297,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 						                  jtom,
 						                  ktom,
 						                  bc_ik|bc_ij);
+
 #ifdef STATIC_LISTS
 					}
 #elif LOWMEM_LISTS
@@ -2973,7 +3100,7 @@ int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	return 0;
 }
 
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 void *visual_atoms(void *ptr) {
 #else
 int visual_atoms(t_moldyn *moldyn) {
@@ -2986,7 +3113,7 @@ int visual_atoms(t_moldyn *moldyn) {
 	t_visual *v;
 	t_atom *atom;
 	t_vb vb;
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 	t_moldyn *moldyn;
 
 	moldyn=ptr;
@@ -3004,7 +3131,7 @@ int visual_atoms(t_moldyn *moldyn) {
 	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 PTHREADS
+#ifndef VISUAL_THREAD
 		return -1;
 #endif
 	}
@@ -3026,7 +3153,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) {
@@ -3072,7 +3201,7 @@ int visual_atoms(t_moldyn *moldyn) {
 
 	close(vb.fd);
 
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 	pthread_exit(NULL);
 
 }