[34] | 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 | } |
---|