1
2 .SUFFIXES:
3 .SUFFIXES: .cpp .o
4
5 #OBJS = energy.o acceleration.o
6
7 .cpp:
8 mpicxx $< -o $@
9
10 .cpp.o:
11 mpicxx -c $<
12
13 all: energy acceleration
14
15 acceleration: acceleration.cpp $(OBJS)
16 mpicxx $@.cpp -o $@
17
18 energy: energy.cpp $(OBJS)
19 mpicxx $@.cpp -o $@
20
21 tar:
22 @tar czf NBabel_PyCpp.tgz Makefile Plummer.in *.h *.cpp *.py
23
24 zip:
25 @zip NBabel.zip Makefile Plummer.in *.h *.cpp *.py
26
27 clean:
28 @$(RM) -f *~ *.o *.pyc energy acceleration NBabel.zip NBabel.tgz
1 #include <iostream>
2 #include <math.h>
3 #include <numeric>
4 #include <cstdlib>
5 using std::atof;
6
7 #include "mpi.h"
8 #include "star.h"
9
10 using namespace std;
11
12 vector<Star> s;
13
14 void calculate_acceleration() {
15
16 for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
17 si->a.assign(3,0);
18 for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
19 vector<real> rij(3);
20 real init = 0.0;
21 for (vector<Star>::iterator sj = si+1; sj != s.end(); ++sj)
22 if(si!=sj) {
23 for (int i = 0; i != 3; ++i)
24 rij[i] = si->r[i]-sj->r[i];
25 real RdotR = inner_product(rij.begin(), rij.end(),
26 rij.begin(), init);
27 real apre = 1./sqrt(RdotR*RdotR*RdotR);
28 for (int i = 0; i != 3; ++i) {
29 si->a[i] -= sj->m*apre*rij[i];
30 sj->a[i] += si->m*apre*rij[i];
31 }
32 }
33 }
34 }
35
36 int main( int argc, char *argv[] )
37 {
38 MPI_Comm parentcomm;
39
40 MPI_Init( &argc, &argv );
41 MPI_Comm_get_parent( &parentcomm );
42 if (parentcomm == MPI_COMM_NULL)
43 {
44 printf("No, I'm not the parent, I am the child.\n");
45 exit(-1);
46 }
47 else
48 {
49
50 MPI::Intercomm parent = MPI::COMM_WORLD.Get_parent();
51
52 int ns;
53 parent.Recv(&ns, 1, MPI::INT, 0, 999);
54 for (int i=0; i<ns; i++)
55 s.push_back(Star());
56
57 int crun;
58 do {
59
60 // printf("I'm the child.\n");
61 parent.Recv(&crun, 1, MPI::INT, 0, 999);
62 // cout << "Child: crun="<<crun<< flush << endl;
63 if (!crun)
64 break;
65
66 // cout << "Receiving " << ns << " stars from parent." << endl << flush ;
67 double scluster[ns*7];
68 // cout << "Child: Receiving stars" << endl;
69 parent.Recv(&scluster, ns*7, MPI::DOUBLE, 0, 999);
70 // cout << "Child: Received stars" << endl;
71
72 for (int i=0; i<ns; i++) {
73 s[i].m = scluster[i*7];
74 for (int j=0; j<3; j++) {
75 s[i].r[j] = scluster[i*7+j+1];
76 s[i].v[j] = scluster[i*7+j+4];
77 }
78 // cout <<"Star="<<s[i]<< endl;
79 }
80
81 calculate_acceleration();
82
83 double acc[3*ns];
84 for (int i=0; i<ns; i++) {
85 //for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
86 for (int j=0; j<3; j++) {
87 acc[i*3+j] = s[i].a[j];
88 }
89 }
90 // cout <<"Child: Sending Stars" << flush << endl;
91 parent.Send(acc, 3*ns, MPI::DOUBLE, 0, 999);
92 // cout <<"Child: Send Stars" << endl;
93
94 }
95 while(1);
96
97 }
98
99 // cout << "Terminating run" << endl;
100 fflush(stdout);
101 MPI_Finalize();
102 return 0;
103 }
1 #!/usr/bin/python
2
3 from __future__ import division
4
5 import operator
6 import pickle
7 from numpy import *
8
9 class Cluster() :
10 def __init__(self, N) :
11 self.m = zeros((N,1))
12 self.r = zeros((N,3))
13 self.v = zeros((N,3))
14 def __iter__(self):
15 return self
16 def __repr__(self) :
17 tt = 'mass:\n%s\n' % (self.m)
18 tt += 'position:\n%s\n' % (self.r)
19 tt += 'velocity:\n%s' % (self.v)
20 return tt
21
22 def energy(self) :
23 KE = 0.5*sum(self.v*self.v*self.m)
24 PE = 0.0
25 for i in range(len(self.r)) :
26 Ri = self.r[i]
27 mi = self.m[i][0]
28 for j in range(i+1, len(self.r)) :
29 Rj = self.r[j]
30 mj = self.m[j][0]
31 Rij = Ri-Rj
32 PE -= mi*mj/sqrt(dot(Rij,Rij))
33 return PE+KE,KE,PE
34
35 if __name__=="__main__":
36
37 fd = open('Plummer.in', 'r')
38 cl = pickle.load(fd)
39 ET,KE,PE = cl.energy()
40 print ET, KE, PE
41
42
1 #include <iostream>
2 #include <math.h>
3 #include <numeric>
4 #include <cstdlib>
5 using std::atof;
6
7 #include "mpi.h"
8 #include "star.h"
9
10 using namespace std;
11
12 vector<Star> s;
13
14 vector<real> energies() {
15 real init = 0;
16 vector<real> E(3), rij(3);
17 E.assign(3,0);
18 for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si)
19 E[1] += 0.5*si->m*inner_product(si->v.begin(), si->v.end(),
20 si->v.begin(), init);
21 E[2] = 0;
22 for (vector<Star>::iterator si = s.begin(); si != s.end(); ++si) {
23 for (vector<Star>::iterator sj = si+1; sj != s.end(); ++sj) {
24 for (int i = 0; i != 3; ++i)
25 rij[i] = si->r[i]-sj->r[i];
26 E[2] -= si->m*sj->m/sqrt(inner_product(rij.begin(), rij.end(),
27 rij.begin(), init));
28 }
29 }
30 E[0] = E[1] + E[2];
31 return E;
32 }
33
34 int main( int argc, char *argv[] )
35 {
36 MPI_Comm parentcomm;
37
38 MPI_Init( &argc, &argv );
39 MPI_Comm_get_parent( &parentcomm );
40 if (parentcomm == MPI_COMM_NULL)
41 {
42 printf("No, I'm not the parent, I am the child.\n");
43 exit(-1);
44 }
45 else
46 {
47
48 MPI::Intercomm parent = MPI::COMM_WORLD.Get_parent();
49
50 int ns;
51 parent.Recv(&ns, 1, MPI::INT, 0, 999);
52 for (int i=0; i<ns; i++)
53 s.push_back(Star());
54
55 int crun;
56 do {
57
58 // printf("I'm the child.\n");
59 parent.Recv(&crun, 1, MPI::INT, 0, 999);
60 // cout << "Child: crun="<<crun<< flush << endl;
61 if (!crun)
62 break;
63
64 // cout << "Receiving " << ns << " stars from parent." << endl << flush ;
65 double scluster[ns*7];
66 // cout << "Child: Receiving stars" << endl;
67 parent.Recv(&scluster, ns*7, MPI::DOUBLE, 0, 999);
68 // cout << "Child: Received stars" << endl;
69
70 for (int i=0; i<ns; i++) {
71 s[i].m = scluster[i*7];
72 for (int j=0; j<3; j++) {
73 s[i].r[j] = scluster[i*7+j+1];
74 s[i].v[j] = scluster[i*7+j+4];
75 }
76 // cout <<"Star="<<s[i]<< endl;
77 }
78
79 vector<real> E(3);
80 E = energies();
81 real energy[3];
82 for (int i=0; i<3; i++)
83 energy[i] = E[i];
84
85 // cout <<"Child: Sending Stars" << flush << endl;
86 parent.Send(energy, 3, MPI::DOUBLE, 0, 999);
87 // cout <<"Child: Send Stars" << endl;
88
89 }
90 while(1);
91
92 }
93
94 // cout << "Terminating run" << endl;
95 fflush(stdout);
96 MPI_Finalize();
97 return 0;
98 }
1 #!/usr/bin/python
2
3 from __future__ import division
4
5 from numpy import *
6 import pickle
7 import Cluster as sc
8
9 import sys
10 from mpi4py import MPI
11 import numpy
12 import time
13
14 def send_cluster(sc, intercomm) :
15
16 nstar = len(sc.m)
17 scluster = numpy.arange(7*nstar, dtype=numpy.double)
18 for i in range(nstar) :
19 scluster[7*i] = sc.m[i]
20 for j in range(3) :
21 scluster[7*i + j+1] = sc.r[i][j]
22 scluster[7*i + j+4] = sc.v[i][j]
23
24 intercomm.Send([scluster,MPI.DOUBLE], dest=0, tag=999)
25
26 def rec_cluster(sc, intercomm) :
27
28 nstar = len(sc.m)
29
30 acc = zeros(shape(sc.r))
31 sacc = numpy.arange(3*nstar, dtype=numpy.double)
32
33 intercomm.Recv([sacc,MPI.DOUBLE], source=0, tag=999)
34
35 for i in range(nstar) :
36 for j in range(3) :
37 acc[i][j] = sacc[i*3 + j]
38
39 return acc
40
41 def rec_energies(sc, intercomm) :
42
43 energies = numpy.arange(3, dtype=numpy.double)
44 intercomm.Recv([energies,MPI.DOUBLE], source=0, tag=999)
45
46 return energies
47
48 def calculate_acceleration(sc, intercomm) :
49
50 crun = numpy.arange(1, dtype=numpy.int32)
51 crun[0] = 1
52 intercomm.Send([crun,MPI.INT], dest=0, tag=999)
53
54 send_cluster(sc, intercomm)
55
56 # Calculate acceleration
57
58 return rec_cluster(sc, intercomm)
59
60 def calculate_energy(sc, intercomm) :
61
62 crun = numpy.arange(1, dtype=numpy.int32)
63 crun[0] = 1
64 intercomm.Send([crun,MPI.INT], dest=0, tag=999)
65
66 send_cluster(sc, intercomm)
67
68 # Calculate energies
69
70 return rec_energies(sc, intercomm)
71
72 if __name__=="__main__":
73
74 np = 1
75 myargv = []
76 for i in range(1, len(sys.argv)):
77 arg = sys.argv[i]
78 if arg == "-n":
79 i = i + 1
80 np = int(sys.argv[i])
81 else:
82 myargv.append(arg)
83
84 fd = open('Plummer.in', 'r')
85 cl = pickle.load(fd)
86
87 comm_energy = MPI.COMM_SELF.Spawn('./energy', args=myargv, maxprocs=np)
88 nstar = numpy.arange(1, dtype=numpy.int32)
89 nstar[0] = len(cl.m)
90 comm_energy.Send([nstar,MPI.INT], dest=0, tag=999)
91
92 ET0,KE0,PE0 = calculate_energy(cl, comm_energy)
93 print "E0 = ", ET0, KE0, PE0
94
95 t = 0.0
96 tend = 1.0
97 dt = 1e-3
98 k = 0
99
100 comm_acc = MPI.COMM_SELF.Spawn('./acceleration', args=myargv, maxprocs=np)
101 comm_acc.Send([nstar,MPI.INT], dest=0, tag=999)
102
103 # Predictor-Corrector Leapfrog
104 acc0 = calculate_acceleration(cl, comm_acc)
105 while t<tend :
106
107 cl.r += dt*cl.v + 0.5*dt*dt*acc0
108 acc1 = calculate_acceleration(cl, comm_acc)
109 cl.v += 0.5*dt*(acc0+acc1)
110 acc0 = acc1
111 t += dt
112 k += 1
113 if k%10==0 :
114 ET,KE,PE = cl.energy()
115 print "t= ", t, " E= ", ET, KE, PE, " dE= ", (ET-ET0)/ET0
116
117 crun = numpy.arange(1, dtype=numpy.int32)
118 crun[0] = 0
119 comm_acc.Send([crun,MPI.INT], dest=0, tag=999)
120 comm_energy.Send([crun,MPI.INT], dest=0, tag=999)
1 #include <iostream>
2 #include <vector>
3
4 typedef double real;
5
6 using namespace std;
7
8 class Star {
9 public:
10 real m;
11 vector<real> r;
12 vector<real> v;
13 vector<real> a;
14 Star() {
15 r.assign(3,0);
16 v.assign(3,0);
17 a.assign(3,0);
18 }
19 Star(real mass, vector<real>pos, vector<real> vel) {
20 m = mass;
21 r = pos;
22 v = vel;
23 }
24 friend ostream & operator << (ostream &so, const Star &si) {
25 so << si.m << " " << si.r[0] << " "<< si.r[1] << " "<< si.r[2] << " "
26 << si.v[0] << " "<< si.v[1] << " "<< si.v[2] << endl;
27 return so;
28 }
29 };