changes in ballistic term & checking for p before getting random value
[physik/nlsop.git] / nlsop_client.c
1 /*
2  * client nlsop code
3  *
4  * author: frank zirkelbach (frank.zirkelbach@physik.uni-augsburg.de)
5  *
6  * this program tries helping to understand the amorphous depuration
7  * and recrystallization of SiCx while ion implantation at temperatures
8  * below 400 degree celsius.
9  * hopefully the program will simulate the stabilization of the
10  * selforganizing lamella structure in the observed behaviour.
11  *
12  * refs: 
13  *  - J. K. N. Lindner. Habil.Schrift, Universitaet Augsburg.
14  *  - Maik Haeberlen. Diplomarbeit, Universitaet Augsburg.
15  *
16  * Copyright (C) 2004 Frank Zirkelbach
17  *
18  * This program is free software; you can redistribute it and/or modify
19  * it under the terms of the GNU General Public License as published by
20  * the Free Software Foundation; either version 2 of the License, or
21  * (at your option) any later version.
22  *
23  * This program is distributed in the hope that it will be useful,
24  * but WITHOUT ANY WARRANTY; without even the implied warranty of
25  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
26  * GNU General Public License for more details.
27  *
28  * You should have received a copy of the GNU General Public License
29  * along with this program; if not, write to the Free Software
30  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
31  *
32  */
33
34 #define _GNU_SOURCE
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <sys/types.h>
39 #include <sys/stat.h>
40 #include <fcntl.h>
41 #include <unistd.h>
42
43 #include <signal.h>
44
45 #include "nlsop.h"
46 #include "dfbapi.h"
47 #include "random.h"
48
49 #include "network.h"
50 #include "event.h"
51
52 #include "nlsop_general.h"
53
54 #define MAKE_AMORPH(N) *(N)|=AMORPH
55 #define MAKE_CRYST(N) *(N)&=~AMORPH
56
57 /* globals */
58
59 char p_file[MAX_CHARS];
60 char n_e_file[MAX_CHARS];
61 char r_file[MAX_CHARS];
62 t_net *gnet;
63 d3_lattice *gd3_l;
64 info *gmy_info;
65 int *gi;
66 unsigned char dc;
67
68 int get_data_and_calc(t_event *event,void *allineed);
69 int nop(t_event *event,void *allineed);
70
71 int usage(char *prog)
72 {
73  puts("usage:");
74  printf("%s -i ip -p port -r/P/n random/profile/neloss file\n",prog);
75  return 1;
76 }
77
78 /*
79  * nlsop internal functions
80  */
81
82 int sputter(d3_lattice *d3_l)
83 {
84  int i,size;
85  int offh,offl;
86
87  size=d3_l->max_x*d3_l->max_y;
88  offl=0;
89  offh=size;
90
91  for(i=0;i<d3_l->max_z-1;i++)
92  {
93   memcpy(d3_l->status+offl,d3_l->status+offh,size);
94   memcpy(d3_l->extra+offl,d3_l->extra+offh,size*sizeof(int));
95   offl=offh;
96   offh+=size;
97  }
98  memset(d3_l->status+offl,0,size);
99  memset(d3_l->extra+offl,0,size);
100
101  return 1;
102 }
103
104 int process_cell(d3_lattice *d3_l,u32 x,u32 y,u32 z,info *my_info,u32 nel_z)
105 {
106  unsigned char *thiz;
107  int *conc;
108  int i,j;
109  int off;
110  double p,q;
111
112  thiz=d3_l->status+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
113  conc=d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
114  //p=my_info->b*nel_z; // energieuebertrag prop zu nukl. bk
115  p=my_info->b*URAND_MAX; // konstanter energieuebertrag
116  for(i=-(my_info->range);i<=my_info->range;i++)
117  {
118   for(j=-(my_info->range);j<=my_info->range;j++)
119   {
120    if(!(i==0 && j==0))
121    {
122     off=((x+d3_l->max_x+i)%d3_l->max_x)+((y+d3_l->max_y+j)%d3_l->max_x)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y;
123     if(*(d3_l->status+off)&AMORPH) p+=my_info->s*(*(d3_l->extra+off))*URAND_MAX/(i*i+j*j);
124    } 
125   }
126  }
127  p+=*conc*my_info->c*URAND_MAX;
128  if(p>=URAND_MAX) MAKE_AMORPH(thiz);
129  else {
130   if(!(*thiz&AMORPH)) {
131    if(get_rand(URAND_MAX)<=p) MAKE_AMORPH(thiz);
132   }
133   else {
134    /* assume 1-p probability */
135    /* also look for neighbours ! */
136    q=(URAND_MAX-p)>0?URAND_MAX-p:0;
137    j=0;
138    j+=(*(d3_l->status+((x+d3_l->max_x+1)%d3_l->max_x)+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0;
139    j+=(*(d3_l->status+((x+d3_l->max_x-1)%d3_l->max_x)+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0;
140    j+=(*(d3_l->status+x+((y+1+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0;
141    j+=(*(d3_l->status+x+((y-1+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0;
142    j+=(*(d3_l->status+x+y*d3_l->max_x+((z+1+d3_l->max_z)%d3_l->max_z)*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0;
143    j+=(*(d3_l->status+x+y*d3_l->max_x+((z-1+d3_l->max_z)%d3_l->max_z)*d3_l->max_x*d3_l->max_y)&AMORPH)?1:0;
144  
145    p+=((q/6)*j);
146    if(p<=URAND_MAX) {
147     if(get_rand(URAND_MAX)>p) MAKE_CRYST(thiz);
148    }
149   }
150  }
151  
152  return 1;
153 }
154
155 int distrib_c(d3_lattice *d3_l,info *my_info,int step,u32 rj_m,u32 *rj_g)
156 {
157  u32 x,y,z;
158  int i,j,k,c;
159  int offset,off;
160  int carry;
161
162  /* put one c ion somewhere in the lattice */
163  x=get_rand(d3_l->max_x);
164  y=get_rand(d3_l->max_y);
165  z=get_rand_reject(d3_l->max_z,rj_m,rj_g);
166  *(d3_l->extra+x+y*d3_l->max_x+z*d3_l->max_x*d3_l->max_y)+=1;
167  (my_info->cc)++;
168
169  if(step%my_info->diff_rate==0)
170  {
171
172  for(i=0;i<d3_l->max_x;i++)
173  {
174   for(j=0;j<d3_l->max_y;j++)
175   {
176    for(k=0;k<d3_l->max_z;k++)
177    {
178     offset=i+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
179     /* case amorph: amorph <- cryst diffusion */
180     if(*(d3_l->status+offset)&AMORPH)
181     {
182      for(c=-1;c<=1;c++)
183      {
184       if(c!=0)
185       {
186        off=((i+d3_l->max_x+c)%d3_l->max_x)+j*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
187        carry=0;
188        if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));
189        if(carry!=0)
190        {
191         *(d3_l->extra+offset)+=carry;
192         *(d3_l->extra+off)-=carry;
193        }
194       }
195      }
196      for(c=-1;c<=1;c++)
197      {
198       if(c!=0)
199       {
200        off=i+((j+c+d3_l->max_y)%d3_l->max_y)*d3_l->max_x+k*d3_l->max_x*d3_l->max_y;
201        carry=0;
202        if(!(*(d3_l->status+off)&AMORPH)) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));
203        if(carry!=0)
204        {
205         *(d3_l->extra+offset)+=carry; 
206         *(d3_l->extra+off)-=carry; 
207        }
208       }
209      }
210      /* diff in z direction */
211      if(k!=0)
212      {
213       off=i+j*d3_l->max_x+(k-1)*d3_l->max_x*d3_l->max_y;
214       carry=0;
215       if(!*(d3_l->status+off)&AMORPH) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));
216       if(carry!=0)
217       {
218        *(d3_l->extra+off)-=carry;
219        *(d3_l->extra+offset)+=carry;
220       }
221      }
222      if(k!=d3_l->max_z-1)
223      {
224       off=i+j*d3_l->max_x+(k+1)*d3_l->max_x*d3_l->max_y;
225       carry=0;
226       if(!*(d3_l->status+off)&AMORPH) carry=(int)(my_info->dr_ac*(*(d3_l->extra+off)));
227       if(carry!=0)
228       {
229        *(d3_l->extra+off)-=carry;
230        *(d3_l->extra+offset)+=carry;
231       }
232      }
233     }
234    } /* for z */
235   } /* for y */
236  } /* for x */
237
238  } /* if step modulo diff_rate == 0 */
239
240  return 1;
241 }
242
243 u32 get_reject_graph(info *my_info,d3_lattice *d3_l,char *file,u32 *graph) {
244  double a,b;
245  int i,j,k;
246  int fd;
247  char buf[32],*p;
248  unsigned char *flag;
249  u32 max;
250
251  max=0;
252  if((fd=open(file,O_RDONLY))<0)
253  {
254   puts("cannot open file to calculate rejection graph");
255   return -1;
256  }
257  if((flag=(unsigned char *)malloc(d3_l->max_z))==NULL)
258  {
259   puts("cannot malloc flag memory for rejection graph");
260   return -1;
261  }
262  memset(flag,0,d3_l->max_z);
263  memset(graph,0,d3_l->max_z*sizeof(u32));
264  /* get fixpoints */
265  k=1;
266  while(k)
267  {
268   for(i=0;i<32;i++)
269   {
270    k=read(fd,&buf[i],1);
271    if((buf[i]=='\n')||(k==0)) break;
272   }
273   if(k)
274   {
275    p=strtok(buf," ");
276    a=atof(p)/10; /* nm */
277    p=strtok(NULL," ");
278    b=atof(p);
279    if(a>d3_l->max_z*CELL_LENGTH) k=0;
280    else 
281    {
282     graph[(int)(a/CELL_LENGTH)]=(int)(URAND_MAX/100*b);
283     flag[(int)(a/CELL_LENGTH)]=1;
284    }
285   }
286  }
287  /* do (linear) interpolation here! */
288  i=0;
289  a=0;
290  while(i<d3_l->max_z)
291  {
292   /* graph[0] is 0! */
293   j=i;
294   i++;
295   while(flag[i]==0&&i<d3_l->max_z) i++;
296   for(k=j+1;k<i;k++) graph[k]=(int)((k-j)*((int)graph[i]-(int)graph[j])/(i-j))+graph[j];
297   if(graph[i]>max) max=graph[i];
298  }
299
300  free(flag);
301
302 #ifdef DEBUG_INTERPOL_PROFILE
303  printf("debug: %s (interpolated profile)\n",file);
304  for(i=0;i<d3_l->max_z;i++) printf("%d %d\n",i,graph[i]);
305 #endif
306
307  return max;
308 }
309
310 void send_data(int signum) {
311
312   int c;
313
314   c=gd3_l->max_x*gd3_l->max_y*gd3_l->max_z;
315
316   network_send_chan(gnet,0,&dc,1);
317   network_send_chan(gnet,0,(unsigned char *)gd3_l,sizeof(d3_lattice));
318   network_send_chan(gnet,0,(unsigned char *)gmy_info,sizeof(info));
319   network_send_chan(gnet,0,gd3_l->status,c*sizeof(unsigned char));
320   network_send_chan(gnet,0,(unsigned char *)gd3_l->extra,c*sizeof(int));
321   network_send_chan(gnet,0,(unsigned char *)gi,sizeof(int));
322 }
323
324
325 /*
326  * main program
327  */
328
329 int main(int argc,char **argv)
330 {
331
332   char server_ip[16];
333   int port;
334   t_net net;
335   t_event event;
336   unsigned char data[256];
337   int i;
338
339   gnet=&net;
340
341   /* default values */
342   strcpy(server_ip,"137.250.82.105");
343   strcpy(p_file,IMP_PROFILE);
344   strcpy(n_e_file,NEL_PROFILE);
345   strcpy(r_file,"");
346   port=1025;
347
348   /* parse/check argv */
349   for(i=1;i<argc;i++) {
350     if(argv[i][0]=='-') {
351       switch(argv[i][1]) {
352         case 'h':
353           usage(argv[0]);
354           return -1;
355         case 'i':
356           strncpy(server_ip,argv[++i],16);
357           break;
358         case 'r':
359           strcpy(r_file,argv[++i]);
360           break;
361         case 'P':
362           strcpy(p_file,argv[++i]);
363           break;
364         case 'n':
365           strcpy(n_e_file,argv[++i]);
366           break;
367         case 'p':
368           port=atoi(argv[++i]);
369           break;
370         default:
371           usage(argv[0]);
372           return -1;
373       }
374     }
375   }
376   if(!strcmp(server_ip,"")) {
377     usage(argv[0]);
378     return -1;
379   }
380
381   /* event init */
382   event_init(&event,1);
383   event_set_timeout(&event,0,0);
384
385   /* connect to server */
386   network_init(&net,1);
387   network_set_connection_info(&net,0,server_ip,port);
388   if(network_connect(&net,0)==N_E_CONNECT) {
389     printf("unable to connect to server, aborting ...\n");
390     return -1;
391   }
392   network_select(&net,0);
393
394   /* tell server: i am a client, i may work for you */
395   data[0]=NLSOP_CLIENT;
396   network_send(net.connection[0].fd,data,1);
397
398   /* wait for job */
399   event_math(net.connection[0].fd,&event,READ,ADD);
400   printf("idle, waiting for jobs ...\n");
401   event_start(&event,NULL,get_data_and_calc,nop);
402
403   return 1;
404 }
405
406 int nop(t_event *event,void *allineed) {
407
408   printf("\ni did a nop :)\n");
409
410   return 1;
411 }
412
413 int get_data_and_calc(t_event *event,void *allineed) {
414
415   d3_lattice d3_l;
416   info my_info;
417   u32 *c_profile;
418   u32 *n_e_loss;
419   u32 ne_max,ip_max;
420   u32 *nel_z;
421   u32 x_c,y_c,z_c;
422   int i,j;
423   int c_step;
424   unsigned char data;
425
426   t_net *net;
427
428   c_step=0;
429   ne_max=0;
430   ip_max=0;
431
432   net=gnet;
433   gd3_l=&d3_l;
434   gmy_info=&my_info;
435   gi=&i;
436   dc=0;
437
438   printf("got a new job ...\n");
439   
440   /* get info (+data) */
441   network_receive(net->connection[0].fd,&data,sizeof(unsigned char));
442   if(data==NLSOP_NJOB || data==NLSOP_CJOB) {
443     network_receive(net->connection[0].fd,(unsigned char *)&d3_l,
444                     sizeof(d3_lattice));
445     network_receive(net->connection[0].fd,(unsigned char *)&my_info,
446                     sizeof(info));
447     c_step=0;
448     j=d3_l.max_x*d3_l.max_y*d3_l.max_z;
449     d3_l.status=(unsigned char *)malloc(j*sizeof(unsigned char));
450     if(d3_l.status==NULL) {
451       printf("status alloc failed\n");
452       return -1;
453     }
454     d3_l.extra=(int *)malloc(j*sizeof(int));
455     if(d3_l.extra==NULL) {
456       printf("extra malloc failed\n");
457       return -1;
458     }
459     if(data==NLSOP_CJOB) {
460       data=DATA_OK;
461       network_receive(net->connection[0].fd,d3_l.status,
462                       j*sizeof(unsigned char));
463       network_send(net->connection[0].fd,&data,sizeof(unsigned char));
464       network_receive(net->connection[0].fd,(unsigned char *)d3_l.extra,
465                       j*sizeof(int));
466       network_send(net->connection[0].fd,&data,sizeof(unsigned char));
467       network_receive(net->connection[0].fd,(unsigned char *)&c_step,
468                       sizeof(int));
469       network_send(net->connection[0].fd,&data,sizeof(unsigned char));
470     }
471   }
472   else {
473     printf("unknown instruction, restarting ...\n");
474     return -1;
475   }
476
477   printf("starting simulation with following parameters:\n");
478   printf("b = %f | c = %f | s = %f\n",my_info.b,my_info.c,my_info.s);
479   printf("diff every %d steps | diff rate = %f\n",my_info.diff_rate,
480                                                   my_info.dr_ac);
481   printf("current step: %d | total steps: %d\n",c_step,my_info.steps);
482   printf("...\n");
483
484   /* care for signals */
485   dc=DC_QUIT;
486   signal(SIGTERM,send_data);
487
488   /* rand init */
489   if(!strcmp(r_file,"")) rand_init(NULL);
490   else rand_init(r_file);
491
492   /* compute graphs for random number rejection method */
493   if((c_profile=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int)))==NULL) {
494     puts("failed allocating memory for carbon profile graph");
495     return -1;
496   }
497   if((n_e_loss=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int)))==NULL) {
498     puts("failed allocating memory for nuclear energy loss graph");
499     return -1;
500   }
501   ip_max=get_reject_graph(&my_info,&d3_l,p_file,c_profile);
502   ne_max=get_reject_graph(&my_info,&d3_l,n_e_file,n_e_loss);
503
504   /* array nel_z[] keeping nuclear energy loss values scaled to URAND_MAX */
505   nel_z=(u32 *)malloc(d3_l.max_z*sizeof(unsigned int));
506   if(nel_z==NULL) {
507     printf("failed allocating nel_z array mem\n");
508     return -1;
509   }
510   for(i=0;i<d3_l.max_z;i++) nel_z[i]=URAND_MAX*(1.0*n_e_loss[i]/ne_max);
511
512   /* this should be obsolete - z is high enough - we check now! */
513   if(c_profile[d3_l.max_z-1]!=0) {
514     printf("max_z (%d) too small - sputtering not possible\n",d3_l.max_z);
515     return -1;
516   }
517
518   /* sputtering really possible ?*/
519   if(n_e_loss[d3_l.max_z-1]!=0)
520     printf("warning: max_z (%d) too small, there may be amorphous volumes\n",
521            d3_l.max_z);
522
523   /*
524    * start of simulation
525    */
526
527   i=(c_step?c_step:0);
528   while(i<my_info.steps) {
529     for(j=0;j<my_info.cpi;j++) {
530       x_c=get_rand(d3_l.max_x);
531       y_c=get_rand(d3_l.max_y);
532       z_c=get_rand_reject(d3_l.max_z,ne_max,n_e_loss);
533       //z_c=get_rand(d3_l.max_z);
534       process_cell(&d3_l,x_c,y_c,z_c,&my_info,nel_z[z_c]);
535     }
536     distrib_c(&d3_l,&my_info,i,ip_max,c_profile);
537     i++;
538     if(i%my_info.save_rate==0) {
539       dc=DC_OK;
540       send_data(0);
541       dc=DC_QUIT;
542     }
543     if(i%my_info.s_rate==0) sputter(&d3_l);
544   }
545
546   /* finished */
547   dc=DC_END;
548   send_data(0);
549   dc=DC_QUIT;
550
551   /* shutdown/free/close everything now ... */
552   free(d3_l.status);
553   free(d3_l.extra);
554   free(c_profile);
555   free(n_e_loss);
556   free(nel_z);
557
558   return 1;
559 }