1 | #include "NBody.hpp" |
---|
2 | |
---|
3 | extern FILE *snapshotsFile, *diagnosticsFile; |
---|
4 | extern void* me; // Pointing to NBody class; |
---|
5 | |
---|
6 | /////////////////////////////////// NBody::put_snapshot() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
---|
7 | // |
---|
8 | // Writes a single snapshot on the output snapshot file. |
---|
9 | // |
---|
10 | /////////////////////////////////// NBody::put_snapshot() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
---|
11 | |
---|
12 | void NBody::put_snapshot() |
---|
13 | { |
---|
14 | if (snapshotsFile == NULL){ |
---|
15 | std::cout << "Error: snapshot file not opened." << std::endl; |
---|
16 | ((NBody*)me)->cleanup(); |
---|
17 | exit(1); |
---|
18 | } |
---|
19 | |
---|
20 | // cout.precision(16); // Full double precision; |
---|
21 | fprintf(snapshotsFile, "%d\n", numParticles); // Total particle number; |
---|
22 | fprintf(snapshotsFile, "%f\n", curr_time_step); // Current time; |
---|
23 | |
---|
24 | for (int i = 0; i < numParticles ; i++){ |
---|
25 | |
---|
26 | fprintf(snapshotsFile, "%f\t", pos[i*4+3]); // Mass of particle i; |
---|
27 | |
---|
28 | for (int k = 0; k < 3; k++) |
---|
29 | fprintf(snapshotsFile, "%f\t", pos[i*4+k]); // Position of particle i; |
---|
30 | |
---|
31 | for (int k = 0; k < 3; k++) |
---|
32 | fprintf(snapshotsFile, "%f\t", vel[i*4+k]); // Velocity of particle i; |
---|
33 | |
---|
34 | fprintf(snapshotsFile, "\n"); |
---|
35 | } |
---|
36 | fprintf(snapshotsFile, "\n"); |
---|
37 | } |
---|
38 | |
---|
39 | |
---|
40 | /////////////////////////////////// NBody::epot_particle() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
---|
41 | // |
---|
42 | // Computes the potential energy of a single particle; |
---|
43 | // |
---|
44 | /////////////////////////////////// NBody::epot_particle() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
---|
45 | |
---|
46 | float NBody::epot_particle(int currPart) |
---|
47 | { |
---|
48 | float ePot = 0.0; |
---|
49 | for (int i=0; i<numParticles; i++){ |
---|
50 | |
---|
51 | if (currPart != i){ |
---|
52 | |
---|
53 | int idx_curr = currPart*4; // The current particle's position index in the "pos" vector; |
---|
54 | int idx_other = i*4; // The other particle's position index in the "pos" vector; |
---|
55 | |
---|
56 | float r[3], r2= 0.0; |
---|
57 | for (int k = 0; k < 3 ; k++){ |
---|
58 | r[k]= pos[idx_other+k] - pos[idx_curr+k]; |
---|
59 | r2+= r[k]*r[k]; |
---|
60 | } |
---|
61 | |
---|
62 | // To avoid Division by Zero: |
---|
63 | if (r2 == 0) r2= 0.1e-50; |
---|
64 | |
---|
65 | // ePot = -M1*M2/r; |
---|
66 | ePot += -pos[idx_curr+3]*pos[idx_other+3]/sqrt(r2); |
---|
67 | } |
---|
68 | } |
---|
69 | |
---|
70 | return ePot; |
---|
71 | } |
---|
72 | |
---|
73 | /////////////////////////////////// NBody::write_diagnostics() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
---|
74 | // |
---|
75 | // Writes diagnostics on the diagnostics file: |
---|
76 | // current time; number of integration steps so far; |
---|
77 | // kinetic, potential, and total energy; absolute and |
---|
78 | // relative energy errors since the start of the run. |
---|
79 | // |
---|
80 | /////////////////////////////////// NBody::write_diagnostics() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
---|
81 | |
---|
82 | void NBody::write_diagnostics() |
---|
83 | { |
---|
84 | |
---|
85 | float Ekin = 0; // The Kinetic energy of the n-body system; |
---|
86 | for (int i = 0; i < numParticles ; i++){ |
---|
87 | |
---|
88 | int idx = i*4; |
---|
89 | |
---|
90 | for (int k = 0; k < 3 ; k++){ |
---|
91 | // pos[idx+3] = mass; |
---|
92 | // Ekin = 1/2*m*(v^2) |
---|
93 | |
---|
94 | Ekin += 0.5 * pos[idx+3] * vel[idx+k] * vel[idx+k]; |
---|
95 | } |
---|
96 | } |
---|
97 | |
---|
98 | float Epot= 0; // The Potential energy of the n-body system; |
---|
99 | for (int i=0; i<numParticles; i++) |
---|
100 | Epot += epot_particle(i); |
---|
101 | |
---|
102 | Epot *= 0.5; // Against double counting; |
---|
103 | |
---|
104 | float Etot = Ekin+Epot; // The total energy of the actual n-body system; |
---|
105 | |
---|
106 | if (initFlag){ // At first pass, the initial energy is equal to the |
---|
107 | Etot_init = Etot; // total energy of the system. We will keep this initial |
---|
108 | // energy for error detections; |
---|
109 | } |
---|
110 | |
---|
111 | if (!initFlag){ // At first pass we don't write any diagnostics; |
---|
112 | if (diagnosticsFile == NULL){ |
---|
113 | std::cout << "Error: diagnostics file not opened." << std::endl; |
---|
114 | ((NBody*)me)->cleanup(); |
---|
115 | exit(1); |
---|
116 | } |
---|
117 | |
---|
118 | // cout.precision(16); // Full double precision; |
---|
119 | fprintf(diagnosticsFile, "At time t = %f, after %ld steps:\n", curr_time_step, taken_steps); |
---|
120 | fprintf(diagnosticsFile, " E_init = %f,\n E_tot = %f,\n E_kin = %f,\n E_pot = %f\n", Etot_init, Etot, Ekin, Epot); |
---|
121 | fprintf(diagnosticsFile, " Absolute energy error: E_tot - E_init = %f\n", Etot - Etot_init); |
---|
122 | fprintf(diagnosticsFile, " Relative energy error: (E_tot - E_init) / E_init = %f\n", (Etot - Etot_init) / Etot_init); |
---|
123 | }else{ |
---|
124 | initFlag = false; |
---|
125 | } |
---|
126 | } |
---|