00001
00002
00003
00004
00005
00006
00007
00008 #include <stdio.h>
00009 #include <sys/types.h>
00010 #ifdef INTERNAL_OUT
00011 # undef INTERNAL_OUT
00012 # define INTERNAL_OUT 1
00013 #else
00014 # define INTERNAL_OUT 0
00015 #endif
00016
00017 unsigned int h3wait_();
00018
00019
00020 #define harpreal double
00021
00022 #define NMAX 1000
00023 #define IPMAX 3
00024 #define JPMAX 26
00025 #define BUFMAX 50
00026 #define NDIM 3
00027
00028 #define NNBMAX 1024
00029
00030 int npipes = IPMAX;
00031 double xjmem[NMAX][NDIM];
00032 harpreal vjmem[NMAX][NDIM];
00033 double xjpred[NMAX][NDIM];
00034 harpreal vjpred[NMAX][NDIM];
00035 harpreal ajmem[NMAX][NDIM];
00036 harpreal jjmem[NMAX][NDIM];
00037 harpreal mjmem[NMAX];
00038 double tjmem[NMAX];
00039
00040 double xjbuf[BUFMAX][JPMAX][NDIM];
00041 harpreal vjbuf[BUFMAX][JPMAX][NDIM];
00042 harpreal ajbuf[BUFMAX][JPMAX][NDIM];
00043 harpreal jjbuf[BUFMAX][JPMAX][NDIM];
00044 harpreal mjbuf[BUFMAX][JPMAX];
00045 double tjbuf[BUFMAX][JPMAX];
00046 int indexbuf[BUFMAX][JPMAX];
00047 int njbuf[BUFMAX];
00048 double timem;
00049
00050 double ximem[IPMAX][NDIM];
00051 harpreal vimem[IPMAX][NDIM];
00052 double hmem[IPMAX];
00053 double epsmem[IPMAX];
00054
00055 int nnb[IPMAX];
00056 int nb[IPMAX][NNBMAX];
00057 static int debug_level = 0;
00058 void h3open_()
00059 {
00060
00061 }
00062
00063 void h3close_()
00064 {
00065 }
00066
00067 int h3npipe_()
00068 {
00069 return npipes;
00070 }
00071
00072 void h3setnboards_(nb)
00073 int * nb;
00074 {
00075
00076 }
00077
00078 int h3getnboards_()
00079 {
00080 return 1;
00081
00082 }
00083
00084 unsigned int h3wait_()
00085 {
00086 return 0;
00087 }
00088
00089 void h3setti_(ti)
00090 double * ti;
00091 {
00092 timem= *ti;
00093 }
00094
00095 int h3jpmax_()
00096 {
00097
00098 return JPMAX;
00099 }
00100
00101
00102 #ifdef G77
00103 #define h3jpdma_indirect_ h3jpdma_indirect__
00104 #define h3mjpdma_indirect_ h3mjpdma_indirect__
00105 #define h3mjpdma_start_ h3mjpdma_start__
00106 #define h3mjpdma_flush_ h3mjpdma_flush__
00107 #endif
00108
00109 void h3setmode_()
00110 {}
00111
00112
00113 void h3jpdma_indirect_(int *,int*,double[][3],double[][3],double[][3],
00114 double[][3],double*,double* ,int*);
00115
00116 void h3mjpdma_indirect_(int *,int*,double[][3],double[][3],double[][3],
00117 double[][3],double*,double* ,int*,int*);
00118 void h3mjpdma_start_(int *);
00119 void h3mjpdma_flush_();
00120
00121 void h3jpdma_indirect_(nj,hostindex,xj,vj,aj,jj,mj,tj,mode)
00122 int * nj;
00123 int hostindex[];
00124 double xj[][3];
00125 double vj[][3];
00126 double aj[][3];
00127 double jj[][3];
00128 double mj[];
00129 double tj[];
00130 int *mode;
00131 {
00132 int zero = 0;
00133 h3mjpdma_indirect_(nj,hostindex,xj,vj,aj,jj,mj,tj,mode,&zero);
00134 h3mjpdma_start_(&zero);
00135 }
00136
00137
00138 static int buff_used = 0;
00139 void h3mjpdma_indirect_(nj,hostindex,xj,vj,aj,jj,mj,tj,mode,buff_id)
00140 int * nj;
00141 int hostindex[];
00142 double xj[][3];
00143 double vj[][3];
00144 double aj[][3];
00145 double jj[][3];
00146 double mj[];
00147 double tj[];
00148 int *mode;
00149 int * buff_id;
00150 {
00151 register int i;
00152 int bufno;
00153 if (debug_level) fprintf(stderr,"(h3mjpdma) called %d %d\n",*nj, *mode);
00154 if(*nj > JPMAX){
00155 fprintf(stderr,"(h3jpdma) Too large nj (%d > %d\n",nj, JPMAX);
00156 exit(-1);
00157 }
00158 if(*buff_id >= 0){
00159 bufno = *buff_id;
00160 }else{
00161 bufno = buff_used;
00162 buff_used ++;
00163 }
00164 njbuf[bufno] = *nj;
00165 for(i=0;i<*nj;i++){
00166 int k;
00167 int hostid;
00168 hostid = hostindex[i]-1;
00169 for(k=0;k<NDIM;k++){
00170 xjbuf[bufno][i][k] = xj[hostid][k];
00171 vjbuf[bufno][i][k] = vj[hostid][k];
00172 ajbuf[bufno][i][k] = aj[hostid][k];
00173 jjbuf[bufno][i][k] = jj[hostid][k];
00174 }
00175 tjbuf[bufno][i] = tj[hostid];
00176 mjbuf[bufno][i] = mj[hostid];
00177 indexbuf[bufno][i] = hostid;
00178 }
00179 }
00180
00181 void h3mjpdma_start_(buff_id)
00182 int * buff_id;
00183 {
00184 register int i,k, j, buff;
00185 if (debug_level> 0) fprintf(stderr,"(h3mjpdma_start) called %d \n",*buff_id);
00186 buff = *buff_id;
00187 for(i=0;i<njbuf[buff];i++){
00188 int k;
00189 int j;
00190 j= indexbuf[buff][i];
00191 if (debug_level > 2){
00192 printf("jpdma_start %d %d %d %f\n",
00193 i, j, buff, xjbuf[buff][i][0]);
00194 }
00195 for(k=0;k<NDIM;k++){
00196 xjmem[j][k] = xjbuf[buff][i][k];
00197 vjmem[j][k] = vjbuf[buff][i][k];
00198 ajmem[j][k] = ajbuf[buff][i][k];
00199 jjmem[j][k] = jjbuf[buff][i][k];
00200 }
00201 tjmem[j]= tjbuf[buff][i];
00202 mjmem[j]= mjbuf[buff][i];
00203 }
00204
00205 }
00206
00207 void h3jpdma_flush_()
00208 {
00209 int i;
00210 for(i=0;i<buff_used; i++){
00211 h3mjpdma_start_(&i);
00212 }
00213 buff_used = 0;
00214 }
00215
00216 void h3jpdma_reset_()
00217 {
00218 buff_used = 0;
00219 }
00220
00221 static int njset;
00222 h3calc_firsthalf_(nj,ni,xi,vi,eps2,h2)
00223 int * nj;
00224 int * ni;
00225 double xi[][3];
00226 double vi[][3];
00227 double eps2[];
00228 double h2[];
00229 {
00230 int i,j,k;
00231 njset = *nj;
00232 if(*ni > IPMAX){
00233 fprintf(stderr,"(h3ipdma) Too large ni (%d > %d\n",ni, IPMAX);
00234 exit(-1);
00235 }
00236 for (i=0;i < *ni; i++){
00237 for(k=0;k<3;k++)ximem[i][k] = xi[i][k];
00238 for(k=0;k<3;k++)vimem[i][k] = vi[i][k];
00239 hmem[i] = 2*h2[i];
00240 epsmem[i] = 2*eps2[i];
00241
00242 }
00243 if (debug_level > 0){
00244 for (i=0; i< *ni; i+= 2){
00245 if ((hmem[i] != hmem[i+1]) ||(hmem[i] != hmem[i+1])){
00246 fprintf(stderr,"h3calc_firsthalf, eps/h2 error for %d\n", i);
00247 }
00248 }
00249 }
00250 for (j=0;j< *nj; j++){
00251 double s, s1, s2;
00252 s = timem - tjmem[j];
00253 s1 = 1.5*s;
00254 s2 = 2.0*s;
00255 for (k=0;k<3;k++){
00256 xjpred[j][k]= ((jjmem[j][k]*s + ajmem[j][k])*s + vjmem[j][k])*s
00257 + xjmem[j][k];
00258 vjpred[j][k]= (jjmem[j][k]*s1 + ajmem[j][k])*s2 + vjmem[j][k];
00259 if (debug_level > 2){
00260 printf("predict %d %d %22.14g %22.14g\n",
00261 j,k,xjpred[j][k], vjpred[j][k]);
00262 }
00263 }
00264 }
00265 }
00266
00267 calculate_force_on_one_particle(i,acc,jerk,pot)
00268 int i;
00269 double acc[][3];
00270 double jerk[][3];
00271 double pot[];
00272 {
00273 register int j, k;
00274 double sqrt();
00275 double xi,yi,zi,vxi,vyi,vzi;
00276 register double dx,dy,dz,dvx,dvy,dvz, dr2, drdv, dr2i, dr3i;
00277 register double ax,ay,az,jx,jy,jz, epsi;
00278 ax=ay=az=jx=jy=jz=0.0;
00279 nnb[i] = 0;
00280 pot[i] = 0.0;
00281 xi=ximem[i][0];
00282 yi=ximem[i][1];
00283 zi=ximem[i][2];
00284 vxi=vimem[i][0];
00285 vyi=vimem[i][1];
00286 vzi=vimem[i][2];
00287 epsi = epsmem[i];
00288 for(j=0;j<njset;j++){
00289 dr2 = epsi;
00290 drdv = 0.0;
00291 dx = xjpred[j][0] - xi;
00292 dy = xjpred[j][1] - yi;
00293 dz = xjpred[j][2] - zi;
00294 dvx = vjpred[j][0] - vxi;
00295 dvy = vjpred[j][1] - vyi;
00296 dvz = vjpred[j][2] - vzi;
00297 dr2 = epsi + dx*dx + dy*dy + dz*dz ;
00298 drdv = dx*dvx + dy*dvy + dz*dvz ;
00299 if (dr2 < hmem[i]){
00300 nb[i][nnb[i]] = j;
00301 nnb[i]++;
00302 }
00303 if (dr2 > 0.0){
00304 dr2i = 1.0/dr2;
00305 dr3i = mjmem[j]*dr2i*sqrt(dr2i);
00306 drdv = 3.*drdv*dr2i;
00307 ax += dx*dr3i;
00308 ay += dy*dr3i;
00309 az += dz*dr3i;
00310 jx += (dvx - dx*drdv)*dr3i;
00311 jy += (dvy - dy*drdv)*dr3i;
00312 jz += (dvz - dz*drdv)*dr3i;
00313 }
00314 }
00315 acc[i][0] = ax;
00316 acc[i][1] = ay;
00317 acc[i][2] = az;
00318 jerk[i][0] = jx;
00319 jerk[i][1] = jy;
00320 jerk[i][2] = jz;
00321 }
00322 h3calc_lasthalf_(ni,acc,jerk,pot)
00323 int * ni;
00324 double acc[][3];
00325 double jerk[][3];
00326 double pot[];
00327 {
00328 register int i, j, k;
00329 if(*ni > IPMAX){
00330 fprintf(stderr,"(h3ipdma) Too large ni (%d > %d\n",ni, IPMAX);
00331 exit(-1);
00332 }
00333 if (debug_level) fprintf(stderr,"lasthalf ni=%d\n",*ni);
00334 for (i=0; i< *ni; i++)calculate_force_on_one_particle(i,acc,jerk,pot);
00335 }
00336 h3calc_lasthalf_old(ni,acc,jerk,pot)
00337 int * ni;
00338 double acc[][3];
00339 double jerk[][3];
00340 double pot[];
00341 {
00342 register int i, j, k;
00343 if(*ni > IPMAX){
00344 fprintf(stderr,"(h3ipdma) Too large ni (%d > %d\n",ni, IPMAX);
00345 exit(-1);
00346 }
00347 if (debug_level) fprintf(stderr,"lasthalf ni=%d\n",*ni);
00348 for (i=0; i< *ni; i++){
00349 double sqrt();
00350 register double dx[3], dv[3], dr2, drdv, dr2i, dr3i;
00351 for(k=0;k<3 ; k++){
00352 acc[i][k] = 0.0;
00353 jerk[i][k] = 0.0;
00354 }
00355 nnb[i] = 0;
00356 pot[i] = 0.0;
00357 for(j=0;j<njset;j++){
00358 dr2 = epsmem[i];
00359 drdv = 0.0;
00360 for(k=0;k<3 ; k++){
00361 dx[k] = xjpred[j][k] - ximem[i][k];
00362 dv[k] = vjpred[j][k] - vimem[i][k];
00363 dr2 += dx[k]*dx[k];
00364 drdv += dx[k]*dv[k];
00365 if (debug_level > 2){
00366 fprintf(stderr,"rijk %d %d %d %g %21.14g %21.14g\n",i, j, k, dx[k],
00367 xjpred[j][k],ximem[i][k]);
00368 }
00369 }
00370 if (dr2 < hmem[i]){
00371 nb[i][nnb[i]] = j;
00372 nnb[i]++;
00373 }
00374 if (debug_level > 1){
00375 fprintf(stderr, "i, j, nnb[i], dr2, h = %d %d %d %g %g\n",
00376 i, j, nnb[i], dr2, hmem[i]);
00377 }
00378 if (dr2 > 0.0){
00379 dr2i = 1.0/dr2;
00380 dr3i = mjmem[j]*dr2i*sqrt(dr2i);
00381 drdv = 3.*drdv*dr2i;
00382 for(k=0;k<3 ; k++){
00383 acc[i][k] += dx[k]*dr3i;
00384 jerk[i][k] += (dv[k] - dx[k]*drdv)*dr3i;
00385 if (debug_level > 1){
00386 fprintf(stderr,"force i,j,k %d %d %d %20.10g \n",i, j, k, acc[i][k]);
00387 }
00388 }
00389 }else{
00390 if (debug_level > 1){
00391 fprintf(stderr,"force skip i,j %d %d \n",i, j);
00392 }
00393 }
00394 }
00395 }
00396 }
00397
00398 h3calc_(nj,ni,xi,vi,eps2,h2, acc, jerk, pot)
00399 int * nj;
00400 int * ni;
00401 double xi[][3];
00402 double vi[][3];
00403 double eps2[];
00404 double h2[];
00405 double acc[][3];
00406 double jerk[][3];
00407 double pot[];
00408
00409 {
00410 h3calc_firsthalf_(nj,ni,xi,vi,eps2,h2);
00411 h3calc_lasthalf_(ni,acc,jerk,pot);
00412 }
00413
00414 void h3nbread_(nboards)
00415 int * nboards;
00416 {
00417
00418 }
00419 #ifdef G77
00420 # define set_debug_level_ set_debug_level__
00421 #endif
00422 void set_debug_level_(int * level)
00423 {
00424 debug_level = *level;
00425 }
00426
00427 void h3setdebuglevel_(int * level)
00428 {
00429 set_debug_level_(level);
00430 }
00431 int compare(i,j)
00432 int *i;
00433 int *j;
00434 {
00435 return *i - *j;
00436 }
00437
00438 int h3nblist_(board, chip, nblist)
00439 int * board;
00440 int * chip;
00441 int * nblist;
00442 {
00443
00444 int i;
00445 for(i=0; i< nnb[*chip]; i++){
00446 nblist[i] = nb[*chip][i];
00447 }
00448 return nnb[*chip];
00449 }
00450
00451
