basic plot

debug/Plot2D_binary.py

@@ -0,0 +1,32 @@
+
+# This program prints Hello, world!
+import sys
+import os
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+
+data_filename=sys.argv[1]
+
+data_filedir='debug/data/'
+data_fileFullpath=data_filedir+data_filename+'.bin'
+fileIO=open(data_fileFullpath,"rb")
+data_All=np.fromfile(fileIO,dtype=np.single,count=-1)
+fileIO.close()
+
+#This is a bit messy but esssentiall the first Nx+1 elements are Nx and x_arr
+# the next Ny+1 element are the Ny and y_arr
+# the rest of the file is matrix to be plotted. 
+# I did it this way so there is only one file that contains everything.
+Nx=int(data_All[0])
+x_vec=data_All[1:Nx+1]
+Ny=int(data_All[Nx+1])
+y_vec=data_All[Nx+2:((Nx+2)+Ny)]
+
+fig=plt.figure()
+ax=fig.add_subplot(1,1,1)
+im=ax.plot(x_vec,y_vec)
+
+filename='debug/plots/' + data_filename +'.pdf'
+plt.savefig(filename)
+print("\n Ploted "+ data_filename+ " using python\n")

src/sph.h

@@ -6,21 +6,16 @@
 #include <vector>
 #include "sphash.hpp"
 
-namespace sph {
-void setup_particles(spatial_hash_map<float, particle> &particles, int n) ;
 
+namespace sph {
+void setup();
 
-float step(spatial_hash_map<float, particle> &particles) ;
-float step(std::vector<particle> &particles);
+void setup_particles(spatial_hash_map<float, particle> &particles, int n) ;
+float step(std::vector<particle> &particles) ;
+float step(sim_data &particles);
 Vec3 grad_pressure(particle i, const std::vector<particle> &locals);
 
-void setup();
-
 
-struct sim_data {
-    std::vector<pair<float>> density_pressure;
-    spatial_hash_map<float, particle> particles;
-};
 
 float W_spiky_kernel(Vec3 xi, Vec3 xj);
 

src/sphash-test.cpp

@@ -113,7 +113,8 @@ TEST_CASE("sphash neighbourhood equivalence") {
 
     for (int i = 0; i < same.size(); i++) {
         auto nb_v = sph::get_local_bf(same, sph::support_radius, i);
-        auto nb_hm = hm.get_locality(same[i].pos);
+        std::vector<particle> nb_hm;
+        hm.get_locality(same[i].pos, nb_hm);
 
         CHECK(nb_v.size() == nb_hm.size());
 
@@ -192,17 +193,19 @@ TEST_CASE("sphash locality-perf") {
     setup_particles(same, sph::h);
     auto hm = setup_hashmap(same);
 
+    std::vector<particle> loc;
+
     {
         Timer t ("sphash-timer");
 
         for (auto p: hm.particles) {
-            hm.get_locality(p.second.pos);
+            hm.get_locality(p.second.pos, loc);
         }
         for (auto p: hm.particles) {
-            hm.get_locality(p.second.pos);
+            hm.get_locality(p.second.pos, loc);
         }
         for (auto p: hm.particles) {
-            hm.get_locality(p.second.pos);
+            hm.get_locality(p.second.pos, loc);
         }
     }
 }

src/sphash.hpp

@@ -2,6 +2,7 @@
 #pragma once
 #include <memory>
 
+#include <cassert>
 #include <cstdlib>
 #include <memory>
 #include <cstddef>
@@ -15,6 +16,9 @@
 
 #include <iostream>
 
+
+struct sim_data;
+
 template <typename V> class aligned_array {
     private:
     std::size_t _size;
@@ -42,10 +46,6 @@ template <typename V> class aligned_array {
     }
 };
 
-template <typename V> class aligned_vector {
-
-};
-
 
 template <typename T, typename V> class spatial_hash_map {
 
@@ -92,10 +92,12 @@ template <typename T, typename V> class spatial_hash_map {
 
         int repair_map() {
             int did = updates.size();
+            int done = 0;
             for (auto p : updates) {
                 for (auto e = particles.begin(); e != particles.end();) {
                     if (e->second.pos == p.second) {
                         e = particles.erase(e);
+                        done++;
                         break;
                     } else {
                         ++e;
@@ -106,13 +108,17 @@ template <typename T, typename V> class spatial_hash_map {
                 insert(p.first.pos, p.first);
             }
             updates.clear();
+            assert(done == did);
             return did;
 
         }
 
-        std::vector<particle> 
-            get_locality(Vec3 pos) 
+        void
+            get_locality(Vec3 pos, std::vector<particle> &loc) 
             {
+                loc.clear();
+
+
             /* Can keep pos key around for deletion hopefully
              * https://en.cppreference.com/w/cpp/container/multimap/find */
 
@@ -162,6 +168,8 @@ template <typename T, typename V> class spatial_hash_map {
                 }
             };
 
+            Vec3 sign = {signs[0], signs[1], signs[2]};
+
             struct mini_vec nums {};
             nums.insert(hash_func(pos));
             for (int i = 0; i < 3; i++) {
@@ -179,9 +187,6 @@ template <typename T, typename V> class spatial_hash_map {
                 }
             }
 
-            std::vector<particle> results;
-            results.reserve(50);
-
             /*
              * For everything in spatially adjacent cells we do the brute force
              * test.
@@ -204,7 +209,7 @@ template <typename T, typename V> class spatial_hash_map {
                                 pos.y - p.pos.y,
                                 pos.z - p.pos.z,
                                 }) <= sph::support_radius) {
-                        results.push_back(p);
+                        loc.push_back(p);
                     }
                 }
 #ifdef debug
@@ -212,8 +217,32 @@ template <typename T, typename V> class spatial_hash_map {
 #endif
         }
 
-        return results;
     }
 };
 
 
+class sim_data {
+    public:
+    spatial_hash_map<float, particle> particles;
+    std::vector<particle_locality> working_area;
+
+    sim_data(float support_radius, int num_particles) :
+        particles{spatial_hash_map<float, particle> {support_radius}} {
+
+        for (int i = 0; i < num_particles; i++) {
+            working_area.push_back(particle_locality {});
+            working_area[i].locality.reserve(AVG_BUCKET_SIZE_EST);
+        }
+
+    }
+
+    void reserve() {
+        for (int i = particles.size(); i < working_area.size(); i++) {
+            working_area.push_back(particle_locality {});
+            working_area[i].locality.reserve(AVG_BUCKET_SIZE_EST);
+        }
+    }
+    sim_data() = delete;
+};
+
+

src/sphfunctions.cpp

@@ -31,6 +31,7 @@ float visc_constant = 0.0007;
 float stiffness_k = 1100;
 float h = 0.1; // smoothing radius
 float mass;
+float gas_const =0.000001; 
 
 float invh;
 float invmass; // 1/mass
@@ -60,6 +61,9 @@ void setup() {
     inv_restdenssity = 1.0/rest_density;
 }
 
+/*
+ * spiky kernel for density , and hence pressure calculation.
+ */
 float spiky_kernel_f(float q) {
     float s = 3.0 / (2 * M_PI);
     if (0 <= q && q < 1) {
@@ -82,6 +86,17 @@ float spiky_kernel_dfdq(float q) {
     }
 }
 
+Vec3 gradW_muller_spiky(Vec3 r) {
+    float q = mod_vec3(r);
+    float res = (-45 * invh * invh * invh*invh * invh * invh / M_PI );
+    if (0 <= q && q <= h) {
+        return r * (1/q) * (h - q) * (h - q) * res;
+    } 
+
+    return Vec3 {0,0,0};
+
+}
+
 float W_spiky_kernel(Vec3 xi, Vec3 xj) {
     float q = mod_vec3(Vec3{
         xi.x - xj.x,
@@ -90,17 +105,52 @@ float W_spiky_kernel(Vec3 xi, Vec3 xj) {
     return spiky_kernel_f(q) * (invh * invh * invh);
 }
 
+float W_muller_poly6_kernel(float r) {
+    float res = 315.0 / 64;
+    res /= M_PI;
+    res *= pow(invh, 9);
+    if (0 <= r && r <= h) {
+        return res * (h * h - r * r) * (h * h - r * r)* (h * h - r * r);
+    }
+    return 0;
+}
+
+Vec3 gradW_muller_poly6(Vec3 r) {
+    float res = 915.0 / (32 * M_PI * pow(h,9));
+    float m = mod_vec3(r);
+
+
+
+}
+
+float W_muller_spiky_kernel(float r) {
+    if (r <= h  && 0 <= r) {
+        return  (15 / M_PI) * invh * invh * invh * invh * invh * invh 
+            * (h - r)* (h - r)* (h - r);
+    }
+    return 0;
+}
+
+
 float W_spiky_kernel(float q) {
     return spiky_kernel_f(q) * (invh * invh * invh);
 }
 
+inline float laplW_viscosity_kernel(float q) {
+
+    if (q <= h && 0 <= q) {
+        return (h - q) *invh * invh * invh * invh * invh * invh * 45 / M_PI;
+    }
+    return 0;
+}
+
 
 float mod_gradW_spiky (Vec3 xi, Vec3 xj) {
     float q = mod_vec3(Vec3 {
         xi.x - xj.x,
         xi.y - xj.y,
-        xi.z - xj.z}) * invh;
-    // TODO: CHECK
+        xi.z - xj.z}) ;
+    // TODO: CHECK / h
     return spiky_kernel_dfdq(q) * (invh * invh * invh * invh);
 }
 
@@ -136,7 +186,7 @@ Vec3 grad_pressure(particle i, const std::vector<particle> &locals) {
     for (const auto &j : locals) {
         float inv_dens_sq = 1.0/(j.density * j.density);
         float temp = mass * (i.pressure * inv_dens_sq  + (j.pressure * inv_dens_sq));
-        Vec3 gw = gradW_spiky(i.pos, j.pos);
+        Vec3 gw = gradW_spiky(i.pos , j.pos);
         res = res + gw * temp;
     }
 
@@ -167,6 +217,57 @@ inline float pressure(float density) {
     //return (rest_density * 45 * 45 / 7) * (pow(density / rest_density, 7) - 1);
 }
 
+inline float gas_pressure(float density) {
+    /* see test case "pow vs multiply" in sphash-test.cpp */
+    return gas_const * (density - rest_density);
+    //return (rest_density * 45 * 45 / 7) * (pow(density / rest_density, 7) - 1);
+}
+
+
+pair<Vec3> 
+muller_update_density_pressure_gradpressure(particle &i, const std::vector<particle> &locals) 
+{
+    // eq 6
+    double density = 0;
+
+    for (const auto &j : locals) {
+        density += mass * W_muller_poly6_kernel(mod_vec3(i.pos- j.pos));
+    }
+
+    i.density = density;
+    i.pressure = pressure(density);
+    Vec3 lapl_vel {0,0,0};
+    Vec3 grad_pres {0,0,0};
+
+    for (int j =0 ;j < locals.size(); j++) {
+        float inv_dens_sq = 1.0/(locals[j].density * locals[j].density);
+        // last frame density/pressure on locals j
+        auto gradW = gradW_muller_spiky(i.pos - locals[j].pos);
+        //float temp = mass * (i.pressure * inv_dens_sq  + (locals[j].pressure * inv_dens_sq));
+        grad_pres = grad_pres + 
+            ((i.pressure + locals[j].pressure) / (2 * locals[j].density)) * gradW;
+
+        Vec3 vij = i.vel - locals[j].vel;
+        Vec3 xij = i.pos - locals[j].pos;
+
+       // lapl_vel = lapl_vel + 2 * (i.vel - locals[j].vel) * 
+       //     (((mass / locals[j].density) * (Vec3_dot(xij, gradW))) / 
+       //      (Vec3_dot(xij, xij) + 0.01 * h * h));
+
+        //lapl_vel = lapl_vel + laplW_viscosity_kernel(mod_vec3(i.pos - locals[j].pos)) 
+         //   * ((locals[j].vel - i.vel) * (1.0/ locals[j].density));
+    }
+
+    grad_pres = 1.0 * grad_pres * mass;
+    lapl_vel = lapl_vel * visc_constant * mass;
+
+    return {grad_pres, lapl_vel};
+
+}
+
+
+
+
 pair<Vec3> 
 update_density_pressure_gradpressure(particle &i, const std::vector<particle> &locals) 
 {
@@ -186,7 +287,7 @@ update_density_pressure_gradpressure(particle &i, const std::vector<particle> &l
     for (int j =0 ;j < locals.size(); j++) {
         float inv_dens_sq = 1.0/(locals[j].density * locals[j].density);
         // last frame density/pressure on locals j
-        auto gradW = gradW_spiky(i.pos, locals[j].pos);
+        auto gradW = gradW_spiky(i.pos , locals[j].pos);
         float temp = mass * (i.pressure * inv_dens_sq  + (locals[j].pressure * inv_dens_sq));
         grad_pres = grad_pres + gradW * temp;
 
@@ -200,8 +301,6 @@ update_density_pressure_gradpressure(particle &i, const std::vector<particle> &l
 
     }
 
-    grad_pres = grad_pres * density;
-
     return {grad_pres, lapl_vel};
 
 }
@@ -309,6 +408,18 @@ Vec3 lapl_vel_w(const particle i, const std::vector<particle> &locals) {
     return result;
 }
 
+Vec3 muller_visc(const particle p, 
+        const std::vector<particle> &locals) {
+        Vec3 Fvisc = {0,0,0};
+        for (auto j : locals) {
+            float q = mod_vec3(p.pos - j.pos);
+            auto lapl_w = laplW_viscosity_kernel(q);
+            Fvisc = Fvisc + ((j.vel - p.vel) * (1.0/j.density)) * lapl_w;
+        }
+        return Fvisc * mass * visc_constant;
+}
+
+
 float step(std::vector<particle> &particles) {
     std::vector<std::vector<particle>> locals;
 
@@ -392,13 +503,22 @@ float step(std::vector<particle> &particles) {
 }
 
 
-float step(spatial_hash_map<float, particle> &particles) {
-    std::vector<std::vector<particle>> locals;
+float step(struct sim_data &data) {
+
+    auto &particles = data.particles;
+
     int rep = particles.repair_map();
 
     // find neighbours
-    for (auto p : particles.particles) {
-        locals.push_back(particles.get_locality(p.second.pos));
+    //
+    {
+    int i = 0; 
+    for (auto &p : particles.particles) {
+        data.working_area[i].p = p.second;
+        data.working_area[i].part_ref = &p.second;
+        particles.get_locality(p.second.pos, data.working_area[i].locality);
+        i++;
+    }
     }
 
     
@@ -419,85 +539,58 @@ float step(spatial_hash_map<float, particle> &particles) {
 
     float dt = TIMESTEP;
 
-    {
-    int i = 0;
-    for (auto &p: particles.particles) {
-        pair<Vec3> gradp_laplvel = sph::update_density_pressure_gradpressure(p.second, locals[i]);
-     //   Vec3 gradp = sph::grad_pressure(p.second, locals[i]);
-        Vec3 Fp = gradp_laplvel.first;
-        Fp.x *= (-mass / p.second.density);
-        Fp.y *= (-mass / p.second.density);
-        Fp.z *= (-mass / p.second.density);
-
-        //Vec3 Fvisc = lapl_vel(p.second, locals[i]);
-        Vec3 Fvisc = gradp_laplvel.second;
-        Fvisc.x *= mass * visc_constant;
-        Fvisc.y *= mass * visc_constant;
-        Fvisc.z *= mass * visc_constant;
+    {int i = 0;
+    for (auto &pp : data.particles.particles) {
+        particle old_p = pp.second;
+        particle new_p = pp.second;
+
+        pair<Vec3> gradp_laplvel =  sph::update_density_pressure_gradpressure(
+                new_p, 
+                data.working_area[i].locality);
+        Vec3 Fp = -mass *  gradp_laplvel.first;
+        Vec3 Fvisc = mass * visc_constant * gradp_laplvel.second;
 
         Vec3 Fother = Vec3{0,0,-9.8 * mass};
 
         float rev = -0.5;
 
+        Vec3 forces = {
+          Fvisc.x + Fp.x + Fother.x,
+          Fvisc.y + Fp.y + Fother.y,
+          Fvisc.z + Fp.z + Fother.z,
+        };
+
+           new_p.vel = new_p.vel + dt * invmass * forces;
+           new_p.pos = new_p.pos + dt * new_p.vel; 
+
         // stop
-        Vec3 newp = p.second.pos;
-        if (newp.z < 0) {
-            newp.z = 0 + (0 - newp.z);
-            p.second.vel.z *= rev;
+        if (new_p.pos.z < 0) {
+            new_p.pos.z = 0 + (0 - new_p.pos.z);
+            new_p.vel.z *= rev;
         }
 
-        if (newp.x < 0) {
-            //newp.x = 0 - eps - particles[i].pos.x;
-            p.second.vel.x *= rev;
+        if (new_p.pos.x < 0) {
+            new_p.vel.x *= rev;
         }
-        if (newp.x > room_width) {
-            // Fother.x -= mass * 1.0 * (-room_width +newp.x) / dt ;
-            //newp.x = room_width + eps - particles[i].pos.x;
-            p.second.vel.x *= rev;
+        if (new_p.pos.x > room_width) {
+            new_p.vel.x *= rev;
         }
-        if (newp.y < 0) {
-            //Fother.y += mass * 1.0 * (0 - newp.y) / dt ;
-            //newp.y = 0 - eps - particles[i].pos.y;
-            p.second.vel.y *= rev;
+        if (new_p.pos.y < 0) {
+            new_p.vel.y *= rev;
         }
-        if (newp.y > room_width) {
-            //newp.y = room_width + eps - particles[i].pos.y;
-            p.second.vel.y *= rev;
+        if (new_p.pos.y > room_width) {
+            new_p.vel.y *= rev;
         }
 
-
-        forces[i] = {
-          Fvisc.x + Fp.x + Fother.x,
-          Fvisc.y + Fp.y + Fother.y,
-          Fvisc.z + Fp.z + Fother.z,
-        };
-
-       p.second.vel = p.second.vel + (forces[i] + forces_old[i]) * invmass * 0.5 * dt;
-       particle newpart {p.second};
-       newpart.pos = newp;
-       newpart.pos = p.second.pos + p.second.vel * dt + (forces[i] * invmass * 0.5 * dt)* dt;
-
-        /*
-        p.second.vel.x += dt * forces[i].x * invmass;
-        p.second.vel.y += dt * forces[i].y * invmass;
-        p.second.vel.z += dt * forces[i].z * invmass;
-
-        particle newpart {p.second};
-        newpart.pos = newp;
-
-        newpart.pos.x += dt * p.second.vel.x;
-        newpart.pos.y += dt * p.second.vel.y;
-        newpart.pos.z += dt * p.second.vel.z;
-        */
-        if (!particles.update(newpart, p.second.pos)) {
-            p.second.pos = newpart.pos;
+        
+        if (!particles.update(new_p, data.working_area[i].part_ref->pos)) {
+             /* If the particle has not crossed a grid boundary we can update it
+             * directly. */
+             *data.working_area[i].part_ref = new_p;
         }
         i++;
     }
     }
-    forces_old = forces;
-
-
     return dt;
 }
 
@@ -531,7 +624,6 @@ void setup_particles(spatial_hash_map<float, particle> &particles, int n) {
             };
             i ++;
         } 
-
     }
 }
 

src/test.cpp

@@ -1,4 +1,5 @@
 
+#include <algorithm>
 #define doplot
 #include <cmath>
 #include <cstdio>
@@ -12,6 +13,7 @@
 #include <thread>
 #include <fstream>
 #include <sstream>
+#include <string.h>
 #include "type.h"
 #include "svector.h"
 #include "sph.h"
@@ -197,6 +199,94 @@ void process_options(int argc, char **argv) {
        }
 }
 
+void SystemCommandCall(std::string command)
+{
+  int commLen = command.length();
+  char *commandchar = new char[commLen + 1];          // declaring character array
+  commandchar = strcpy(commandchar, command.c_str()); // copying the contents of the string to char array
+  system(commandchar);                                // Creates the directory incase it does not exist
+  delete[] commandchar;
+}
+
+
+void BinarySaveAndPlot2D(int Narr, float *arrX, float *arrY, std::string OutputFileName)
+{
+  float Nf = (float)Narr;
+  
+  //Make a file to write to and open it 
+  std::fstream iofile; // declare file pointer
+  iofile.open("debug/data/" + OutputFileName + ".bin", std::ios::out | std::ios_base::binary);//Open file
+
+  //Write data to file
+  iofile.write((char *)&Nf, sizeof(float));
+  iofile.write((char *)arrX, sizeof(float) * Narr);
+  iofile.write((char *)&Nf, sizeof(float));
+  iofile.write((char *)arrY, sizeof(float) * Narr);
+  iofile.close(); // This is really import as we are going to access the file in python and it needs to closed for python to access it.
+
+  //Call Python to plot data 
+  std::string command = "python debug/Plot2D_binary.py " + OutputFileName;
+  SystemCommandCall(command);
+}
+
+
+struct plotdata_collec {
+    float* time;
+    float* sim_time;
+    float* kinetic_energy;
+    float* potential_energy;
+    float* max_density;
+    float* avg_density;
+    size_t nmemb = 0;
+
+    void reserve(size_t n) {
+        time = (float*)calloc(n, sizeof(float));
+        sim_time = (float*)calloc(n, sizeof(float));
+        kinetic_energy = (float*)calloc(n, sizeof(float));
+        potential_energy = (float*)calloc(n, sizeof(float));
+        max_density = (float*)calloc(n, sizeof(float));
+        avg_density = (float*)calloc(n, sizeof(float));
+    }
+};     
+
+void 
+measure_particle_properties(sim_data &data, float *kinetic, 
+        float *gpotential, float *max_density, float *avg_density) 
+{
+    float kin = 0;
+    float pot = 0;
+    float max_dens = 0;
+    float avg_dens = 0;
+
+    for (auto p: data.particles.particles) {
+        kin +=  Vec3_dot(p.second.vel, p.second.vel) * 0.5  * sph::mass;
+        pot += sph::mass * 9.81 * p.second.pos.z;
+        avg_dens += p.second.density;
+        if (p.second.density > max_dens) 
+            max_dens = p.second.density;
+    }
+
+    *kinetic = kin;
+    *gpotential = pot;
+    *max_density = max_dens;
+    *avg_density = avg_dens / data.particles.size();
+}
+
+void
+add_data_frame(sim_data &data, plotdata_collec &plot, float sim_now, float real_now) 
+{
+    plot.time[plot.nmemb] = sim_now;
+    plot.sim_time[plot.nmemb] = real_now;
+
+    measure_particle_properties(data, 
+            &plot.kinetic_energy[plot.nmemb], 
+            &plot.potential_energy[plot.nmemb],
+            &plot.max_density[plot.nmemb],
+            &plot.avg_density[plot.nmemb]);
+
+    plot.nmemb++;
+}
+
 
 int main(int argc, char **argv) {
     process_options(argc, argv);
@@ -208,13 +298,23 @@ int main(int argc, char **argv) {
     bool once = false;
     float sim_time = 0;
     sph::setup();
-    spatial_hash_map<float, particle> particles {sph::support_radius};
-    sph::setup_particles(particles, options.num_particles);
     auto start = chrono::steady_clock::now();
 
+    sim_data data {sph::support_radius, options.num_particles};
+
+    sph::setup_particles(data.particles, options.num_particles);
+    data.reserve();
+
+    for (auto p: data.particles.particles) {
+        data.working_area.push_back({p.second});
+    }
+
+    plotdata_collec plot_data;
+    plot_data.reserve(options.sim_for / sph::TIMESTEP);
+
 #ifdef doplot
     auto animationFigure = matplot::figure();
-    mkplot(animationFigure, particles);
+    mkplot(animationFigure, data.particles);
     matplot::xlim({-0.1,sph::room_width + 0.1});
     matplot::ylim({-0.1,sph::room_width + 0.1});
     matplot::zlim({-0.1,2});
@@ -222,24 +322,38 @@ int main(int argc, char **argv) {
 
     Timer t {"main"};
 
+    auto last = std::chrono::high_resolution_clock::now();
+    auto now = last;
 
     {
         Timer t {"simulation time"};
     while (sim_time < options.sim_for) {
-        dt = sph::step(particles);
+        dt = sph::step(data);
         frame += dt;
         sim_time += dt;
+
+        now = std::chrono::high_resolution_clock::now();
+        auto delta = now - last;
+        last = now;
+        add_data_frame(data, plot_data, sim_time, delta.count());
+
         if (frame > options.frame_step) {
             std::cout << "Wrote "  << frame << " now " << sim_time << std::endl;
-            write_frame(options.save_path, frames++, frame, particles);
+            write_frame(options.save_path, frames++, frame, data.particles);
 #ifdef doplot
-            mkplot(animationFigure, particles);
+            mkplot(animationFigure, data.particles);
 #endif
             frame = 0;
         }
     }
      }
-    print_lost(particles);
+    print_lost(data.particles);
+
+    BinarySaveAndPlot2D(plot_data.nmemb, plot_data.time, plot_data.max_density, "max_density");
+    BinarySaveAndPlot2D(plot_data.nmemb, plot_data.time, plot_data.avg_density, "avg_density");
+    BinarySaveAndPlot2D(plot_data.nmemb, plot_data.time, plot_data.sim_time, "sim_time");
+    BinarySaveAndPlot2D(plot_data.nmemb, plot_data.time, plot_data.kinetic_energy, "kinetic");
+    BinarySaveAndPlot2D(plot_data.nmemb, plot_data.time, plot_data.potential_energy, "potential");
 
     return 0;
 }

src/type.h

@@ -1,6 +1,7 @@
 
 #pragma once
 
+#define AVG_BUCKET_SIZE_EST 50
 #define ALIGN 1024
 #include <vector>
 
@@ -18,6 +19,13 @@ struct pair {
     T second;
 };
 
+struct particle_locality {
+    particle p;
+    particle *part_ref;
+    std::vector<particle> locality;
+};
+
+
 bool operator==(const particle a, const particle b);
 
 namespace sph {