00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_interaction.h"
00026 #include "orsa_secure_math.h"
00027 #include "orsa_universe.h"
00028
00029 #include <cmath>
00030 #include <iostream>
00031
00032 #include <mpi.h>
00033
00034 using namespace std;
00035
00036 namespace orsa {
00037
00038
00039
00040 Newton_MPI::Newton_MPI() : Interaction() {
00041 av = av_local = bv = 0;
00042 }
00043
00044 Newton_MPI::Newton_MPI(const Newton_MPI &) : Interaction() {
00045 av = av_local = bv = 0;
00046 }
00047
00048 Newton_MPI::~Newton_MPI() {
00049 free(bv);
00050 free(av);
00051 free(av_local);
00052 }
00053
00054 Interaction * Newton_MPI::clone() const {
00055 return new Newton_MPI(*this);
00056 }
00057
00058 void Newton_MPI::Acceleration(const Frame &fr, vector<Vector> &a) {
00059
00060 int rank, size;
00061 MPI_Comm_rank(MPI_COMM_WORLD, &rank );
00062 MPI_Comm_size(MPI_COMM_WORLD, &size );
00063
00064
00065
00066 unsigned int i,j;
00067
00068 Vector d;
00069
00070 double ls;
00071
00072
00073 MPI_Barrier(MPI_COMM_WORLD);
00074
00075 unsigned int frame_size=0;
00076 if (rank==0) {
00077 frame_size = fr.size();
00078 }
00079
00080 MPI_Bcast(&frame_size, 1, MPI_INT, 0, MPI_COMM_WORLD);
00081
00082 if (rank==0) {
00083 g = GetG();
00084 }
00085
00086 MPI_Bcast(&g, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00087
00088
00089
00090 if (frame_size < 2) return;
00091
00092 a.resize(frame_size);
00093 for (i=0;i<frame_size;++i)
00094 a[i].Set(0.0,0.0,0.0);
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 bv = (double *)realloc(bv,4*frame_size*sizeof(double));
00109
00110 if (rank==0) {
00111 unsigned int k=0;
00112 while (k<frame_size) {
00113 bv[4*k] = fr[k].mass();
00114
00115 bv[4*k+1] = fr[k].position().x;
00116 bv[4*k+2] = fr[k].position().y;
00117 bv[4*k+3] = fr[k].position().z;
00118
00119 ++k;
00120 }
00121 }
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 MPI_Bcast(bv, 4*frame_size, MPI_DOUBLE, 0, MPI_COMM_WORLD);
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153 for (i=1+rank;i<frame_size;i+=size) {
00154
00155 for (j=0;j<i;++j) {
00156
00157
00158
00159
00160 if ((bv[4*i]==0.0) && (bv[4*j]==0.0)) continue;
00161
00162
00163 d.Set(bv[4*j+1]-bv[4*i+1],
00164 bv[4*j+2]-bv[4*i+2],
00165 bv[4*j+3]-bv[4*i+3]);
00166
00167
00168
00169 ls = d.LengthSquared();
00170
00171
00172
00173 if (d.IsZero()) {
00174 cerr << "*** Warning: two objects in the same position!" << endl;
00175 continue;
00176 }
00177
00178 d *= secure_pow(ls,-1.5);
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 a[i] = a[i] + d * bv[4*j];
00189 a[j] = a[j] - d * bv[4*i];
00190
00191
00192 }
00193 }
00194
00195
00196
00197
00198 av = (double *)realloc(av,3*frame_size*sizeof(double));
00199
00200 if (rank == 0) {
00201 unsigned int k=0;
00202 while (k<3*frame_size) {
00203 av[k] = 0.0;
00204 ++k;
00205 }
00206 }
00207
00208
00209 av_local = (double *)realloc(av_local,3*frame_size*sizeof(double));
00210
00211 {
00212 unsigned int k=0;
00213 while (k<frame_size) {
00214 av_local[3*k] = a[k].x;
00215 av_local[3*k+1] = a[k].y;
00216 av_local[3*k+2] = a[k].z;
00217 ++k;
00218 }
00219 }
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 MPI_Reduce(av_local, av, 3*frame_size, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00232
00233 if (rank == 0) {
00234 unsigned int k=0;
00235 while (k<frame_size) {
00236 a[k].x = av[3*k];
00237 a[k].y = av[3*k+1];
00238 a[k].z = av[3*k+2];
00239
00240 a[k] *= g;
00241
00242
00243
00244 ++k;
00245 }
00246 }
00247 }
00248
00249 double Newton_MPI::PotentialEnergy(const Frame &f) {
00250 Newton newton;
00251 return (newton.PotentialEnergy(f));
00252 }
00253
00254 }