/******************************************************** % % Written by: % -- % John L. Weatherwax 2006-08-06 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- Compute the short range forces Short range forces include: i) gravitational ii) lenard-jones ... NOT implemented ... */ #include #include #include "mpi.h" #include "global_grid.h" #include "sim_consts.h" #include "computeShortRangeForces.h" /*#define _CONST_FORCE_ 1*/ #define _REAL_FORCE_ 1 /* the gravitational constant */ /*#define gravG 6.67e-11*/ #define gravG 1.0 /* floor on the particle distance */ double DISTANCE_LOWER_BOUND=0.01; /* function declarations: */ int inHomeBox_ONE_NBR(double,double,int); int inHomeBox_TWO_NBR(double,double,int,int); /* creates arrays for data that will be sent to each processor */ void computeShortRangeForces(int tsi){ int ii, cn, cno, jj; double dist, f12, r12_x, r12_y, px, py, pxo, pyo, mx, my, dx, dy; FILE *Fp; /* open a file to save all force calculations into ... */ sprintf(procInfo.pFilename,"ParSimResults/force_calc_proc_%d_ts_%d.dat",procInfo.my_rank,tsi); Fp=fopen(procInfo.pFilename,"w"); /* Compute the short range forces for my original "home" particles. */ for( ii = 0; ii < procInfo.NbPP; ii++ ){ for( jj = ii + 1; jj < procInfo.NbPP; jj++ ){ /* The vector from body "ii" to body "jj" */ dx = procInfo.px[jj] - procInfo.px[ii]; dy = procInfo.py[jj] - procInfo.py[ii]; /*dist = sqrt( dx*dx + dy*dy );*/ dist = sqrtf( dx*dx + dy*dy ); /* the particles are too far apart ... no force update */ if( dist >= simConsts.R ) continue; if (dist < 0.01) { /*printf("particles too close together...mollifying\n"); */ dist += DISTANCE_LOWER_BOUND; } /* The magnitude of the force ... assuming all particles have unit mass */ /* f12 = gravG * procInfo.m[ii] * procInfo.m[jj] / (dist*dist); */ f12 = gravG / (dist * dist); /* The {\em unit} vector from body "ii" to body "jj" */ r12_x = dx / dist; r12_y = dy / dist; #ifdef _CONST_FORCE_ procInfo.fx[ii] += 1.0/3.0; procInfo.fy[ii] += 1.0/5.0; procInfo.fx[jj] += 1.0/3.0; procInfo.fy[jj] += 1.0/5.0; #endif #ifdef _REAL_FORCE_ procInfo.fx[ii] += f12 * r12_x; /* a direct update ... the force on body i due to body j */ procInfo.fy[ii] += f12 * r12_y; procInfo.fx[jj] += -f12 * r12_x; /* a symmetric update ... the force on body j due to body i */ procInfo.fy[jj] += -f12 * r12_y; #endif /* save intermediates ... */ if( procInfo.px[ii] < procInfo.px[jj] ){ /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", procInfo.px[ii], procInfo.py[ii], procInfo.px[jj], procInfo.py[jj], dx,dy,dist,f12,r12_x,r12_y,f12*r12_x,f12*r12_y);*/ fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", procInfo.px[ii], procInfo.py[ii], procInfo.px[jj], procInfo.py[jj], f12*r12_x,f12*r12_y); }else{ /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", procInfo.px[jj], procInfo.py[jj], procInfo.px[ii], procInfo.py[ii], -dx,-dy,dist,f12,-r12_x,-r12_y,-f12*r12_x,-f12*r12_y);*/ fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", procInfo.px[jj], procInfo.py[jj], procInfo.px[ii], procInfo.py[ii], -f12*r12_x,-f12*r12_y); } } } /* some defensive programming */ ii=-1; jj=-1; cn=-1; /* Compute the short range forces between the imported particles and the home particles */ for( ii = 0; ii < procInfo.NbPP; ii++ ){ /* for each particle in the "home" box ... */ for( cn = 0; cn < 9; cn++ ){ /* for each of my compass neighbors I RECIEVED particles FROM (except myself) ... */ if (cn == 4) continue; for (jj = 0; jj < procInfo.numToReceive[cn]; jj++) { /* for each imported particle from this compass neighbor ... */ px = procInfo.partDataToReceive[cn][2 * jj ]; py = procInfo.partDataToReceive[cn][2 * jj + 1]; /* compute the midpoint between these two particles */ mx = 0.5*(procInfo.px[ii] + px); my = 0.5*(procInfo.py[ii] + py); /* this box does not contain the midpoint ... to avoid double counting we don't perform a force update */ if( !inHomeBox_ONE_NBR(mx,my,cn) ) continue; /* The vector from body "ii" to body "jj" */ dx = px - procInfo.px[ii]; dy = py - procInfo.py[ii]; /*dist = sqrt( dx*dx + dy*dy );*/ dist = sqrtf( dx*dx + dy*dy ); /* the particles are too far apart ... no force update */ if( dist >= simConsts.R ) continue; if (dist < 0.01) { /*printf("particles too close together...mollifying\n"); */ dist += DISTANCE_LOWER_BOUND; } /* The magnitude of the force ... assuming all particles have unit mass */ f12 = gravG / (dist * dist); /* The {\em unit} vector from body "ii" to body "p" ("jj") */ r12_x = dx / dist; r12_y = dy / dist; #ifdef _CONST_FORCE_ procInfo.fx[ii] += 1.0/3.0; procInfo.fy[ii] += 1.0/5.0; procInfo.forcesOnReceivedPart[cn][2 * jj ] += 1.0/3.0; procInfo.forcesOnReceivedPart[cn][2 * jj + 1] += 1.0/5.0; #endif #ifdef _REAL_FORCE_ procInfo.fx[ii] += f12 * r12_x; /* a direct update */ procInfo.fy[ii] += f12 * r12_y; procInfo.forcesOnReceivedPart[cn][2 * jj ] += -f12 * r12_x; /* a symmetric update */ procInfo.forcesOnReceivedPart[cn][2 * jj + 1] += -f12 * r12_y; #endif /* save intermediates ... smallest x coordinate ordered first */ if( procInfo.px[ii] < px ){ /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", procInfo.px[ii], procInfo.py[ii], px,py, dx,dy,dist,f12,r12_x,r12_y,f12*r12_x,f12*r12_y);*/ fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", procInfo.px[ii], procInfo.py[ii], px,py, f12*r12_x,f12*r12_y); }else{ /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", px,py, procInfo.px[ii], procInfo.py[ii], -dx,-dy,dist,f12,-r12_x,-r12_y,-f12*r12_x,-f12*r12_y);*/ fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", px,py, procInfo.px[ii], procInfo.py[ii], -f12*r12_x,-f12*r12_y); } } /* for each imported particle ... */ } /* for each compass neighbor except myself ... */ } /* for each particle in the "home" box ... */ /* some defensive programming */ ii=-1; jj=-1; cn=-1; /* Compute the short range forces between two sets of imported particles */ for( cn = 0; cn < 9; cn++ ){ /* for each of my compass neighbors except myself ... */ if( cn == 4 ) continue; for( cno = cn+1; cno < 9; cno++ ){ /* for each of my compass neighbors except myself ... */ /*for( cno = 0; cno < 9; cno++ ){*/ if( cno == 4 ) continue; for( ii = 0; ii < procInfo.numToReceive[cn]; ii++ ){ /* for each imported particle ... */ px = procInfo.partDataToReceive[cn][2 * ii ]; py = procInfo.partDataToReceive[cn][2 * ii + 1]; for( jj = 0; jj < procInfo.numToReceive[cno]; jj++ ){ /* for each imported particle ... */ pxo = procInfo.partDataToReceive[cno][2 * jj ]; pyo = procInfo.partDataToReceive[cno][2 * jj + 1]; /* compute the midpoint between these two particles */ mx = 0.5*(pxo + px); my = 0.5*(pyo + py); /* this box does not contain the midpoint ... to avoid double counting we don't perform a force update */ if( !inHomeBox_TWO_NBR(mx,my,cn,cno) ) continue; /* The vector from body "ii" to body "jj" */ dx = pxo - px; dy = pyo - py; /*dist = sqrt( dx*dx + dy*dy );*/ dist = sqrtf( dx*dx + dy*dy ); /* the particles are too far apart ... no force update */ if( dist >= simConsts.R ) continue; if( dist < 0.01 ){ /*printf("particles too close together...mollifying\n"); */ dist += DISTANCE_LOWER_BOUND; } /* The magnitude of the force ... assuming all particles have unit mass */ f12 = gravG / (dist * dist); /* The {\em unit} vector from body "p"=="ii" to body "po"=="jj" */ r12_x = dx / dist; r12_y = dy / dist; #ifdef _CONST_FORCE_ procInfo.forcesOnReceivedPart[cn][2 * ii ] += 1.0/3.0; procInfo.forcesOnReceivedPart[cn][2 * ii + 1] += 1.0/5.0; procInfo.forcesOnReceivedPart[cno][2 * jj ] += 1.0/3.0; procInfo.forcesOnReceivedPart[cno][2 * jj + 1] += 1.0/5.0; #endif #ifdef _REAL_FORCE_ procInfo.forcesOnReceivedPart[cn][2 * ii ] += f12 * r12_x; /* a direct update */ procInfo.forcesOnReceivedPart[cn][2 * ii + 1] += f12 * r12_y; procInfo.forcesOnReceivedPart[cno][2 * jj ] += -f12 * r12_x; /* a symmetric update */ procInfo.forcesOnReceivedPart[cno][2 * jj + 1] += -f12 * r12_y; #endif /* save intermediates ... */ if( px < pxo ){ /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", px,py pxo,pyo, dx,dy,dist,f12,r12_x,r12_y,f12*r12_x,f12*r12_y);*/ fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", px,py, pxo,pyo, f12*r12_x,f12*r12_y); }else{ /*fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", pxo,pyo, px,py, -dx,-dy,dist,f12,-r12_x,-r12_y,-f12*r12_x,-f12*r12_y);*/ fprintf(Fp,"%20.8f %20.8f %20.8f %20.8f %20.8f %20.8f\n", pxo,pyo, px,py, -f12*r12_x,-f12*r12_y); } } /* for each imported particle ... */ } /* for each imported particle ... */ } /* for each compass neighbor except myself ... */ } /* for each compass neighbor except myself ... */ fclose(Fp); } /* compute whether or not I handle the force update when ONE particles comes from my neighbors */ int inHomeBox_ONE_NBR(double mx,double my,int cn){ int in_x_dom, in_y_dom; in_x_dom = ( procInfo.gdSW[0] < mx ) && ( mx <= procInfo.gdSE[0] ); in_y_dom = ( procInfo.gdSW[1] <= my ) && ( my < procInfo.gdNW[1] ); if( in_x_dom && in_y_dom ){ return 1; /* I should compute this force */ }else{ return 0; /* my neighbor should compute this force */ } } /* compute whether or not I handle the force update when BOTH particles come from my neighbors */ int inHomeBox_TWO_NBR(double mx,double my,int cn,int cno){ int in_x_dom, in_y_dom; in_x_dom = ( procInfo.gdSW[0] < mx ) && ( mx <= procInfo.gdSE[0] ); in_y_dom = ( procInfo.gdSW[1] <= my ) && ( my < procInfo.gdNW[1] ); if( in_x_dom && in_y_dom ){ return 1; /* I should compute this force */ }else{ return 0; /* my neighbor should compute this force */ } }