Templatized SPHSmooth. Reduced memory consumption of simple3DFilter

sample/simple3DFilter.cpp

@@ -7,6 +7,8 @@
 #include <H5Cpp.h>
 #include "hdf5_array.hpp"
 #include <iostream>
+#include <boost/format.hpp>
+#include <boost/bind.hpp>
 
 using namespace std;
 using namespace CosmoTool;
@@ -17,8 +19,13 @@ struct VCoord{
   float v[3];
 };
 
-template<int i>
+using boost::format;
-ComputePrecision getVelocity(const VCoord& v)
+using boost::str;
+typedef boost::multi_array<float, 2> array_type;
+typedef boost::multi_array<float, 3> array3_type;
+typedef boost::multi_array<float, 4> array4_type;
+
+ComputePrecision getVelocity(const VCoord& v, int i)
 {
   return v.v[i];
 }
@@ -32,11 +39,51 @@ typedef SPHSmooth<VCoord> MySmooth;
 typedef MySmooth::SPHTree MyTree;
 typedef MyTree::Cell MyCell;
 
+template<typename FuncT>
+void computeInterpolatedField(MyTree *tree1, double boxsize, int Nres, double cx, double cy, double cz,
+                              array3_type& bins, array3_type& arr, FuncT func, double rLimit2)
+{
+#pragma omp parallel
+  {
+    MySmooth smooth1(tree1, N_SPH);
+    
+#pragma omp for schedule(dynamic) 
+    for (int rz = 0; rz < Nres; rz++)
+      {
+        double pz = (rz)*boxsize/Nres-cz;
+
+        cout << format("[%d] %d / %d") % omp_get_thread_num() % rz % Nres << endl;
+        for (int ry = 0; ry < Nres; ry++)
+          {
+            double py = (ry)*boxsize/Nres-cy;
+            for (int rx = 0; rx < Nres; rx++)
+              {
+                double px = (rx)*boxsize/Nres-cx;
+                
+                MyTree::coords c = { px, py, pz };
+
+                double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+                if (r2 > rLimit2)
+                  {
+                    arr[rx][ry][rz] = 0;
+                    continue;
+                  }
+
+                uint32_t numInCell = bins[rx][ry][rz];
+                if (numInCell > N_SPH)
+                  smooth1.fetchNeighbours(c, numInCell);
+                else
+                  smooth1.fetchNeighbours(c);
+
+                arr[rx][ry][rz] = smooth1.computeSmoothedValue(c, func);
+              }
+          }
+      }
+  }
+}
+
 int main(int argc, char **argv)
 {
-  typedef boost::multi_array<float, 2> array_type;
-  typedef boost::multi_array<float, 3> array3_type;
-  typedef boost::multi_array<float, 4> array4_type;
 
   char *fname1, *fname2;
   double rLimit, boxsize, rLimit2, cx, cy, cz;
@@ -105,8 +152,6 @@ int main(int argc, char **argv)
 
   cout << "Creating smoothing filter..." << endl;
 
-  array3_type out_den_1(boost::extents[Nres][Nres][Nres]);
-  array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]);
 //  array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
  
   cout << "Weighing..." << endl;
@@ -149,55 +194,18 @@ int main(int argc, char **argv)
   }
   
   cout << "Interpolating..." << endl;
-#pragma omp parallel
-  {
-    MySmooth smooth1(&tree1, N_SPH);
 
-#pragma omp for schedule(dynamic) 
+  array3_type interpolated(boost::extents[Nres][Nres][Nres]);
-    for (int rz = 0; rz < Nres; rz++)
-      {
-        double pz = (rz)*boxsize/Nres-cz;
   
-        cout << rz << " / " << Nres << endl;
+  computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz,
-        for (int ry = 0; ry < Nres; ry++)
+                           bins, interpolated, getUnity, rLimit2);
-          {
+  hdf5_write_array(out_f, "density", interpolated);
-            double py = (ry)*boxsize/Nres-cy;
+  //out_f.flush();
-            for (int rx = 0; rx < Nres; rx++)
+  for (int i = 0; i < 3; i++) {
-              {
+      computeInterpolatedField(&tree1, boxsize, Nres, cx, cy, cz,
-                double px = (rx)*boxsize/Nres-cx;
+                               bins, interpolated, boost::bind(getVelocity, _1, i), rLimit2);
-                
+      hdf5_write_array(out_f, str(format("p%d") % i), interpolated);
-                MyTree::coords c = { px, py, pz };
-
-                double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
-                if (r2 > rLimit2)
-                  {
-                    out_v3d_1[rx][ry][rz][0] = 0;
-                    out_v3d_1[rx][ry][rz][1] = 0;
-                    out_v3d_1[rx][ry][rz][2] = 0;
-                    out_den_1[rx][ry][rz] = 0;
-                    //out_rad_1[rx][ry][rz] = -1;
-                    continue;
   }
   
-                uint32_t numInCell = bins[rx][ry][rz];
-                if (numInCell > N_SPH)
-                  smooth1.fetchNeighbours(c, numInCell);
-                else
-                  smooth1.fetchNeighbours(c);
-
-                //out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen();
-                out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>);
-                out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>);
-                out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>);
-                out_den_1[rx][ry][rz] = smooth1.computeSmoothedValue(c, getUnity);
-              }
-          }
-      }
-  }
-      
-  //hdf5_write_array(out_f, "radii", out_rad_1);
-  hdf5_write_array(out_f, "velocity", out_v3d_1);
-  hdf5_write_array(out_f, "density", out_den_1);
-
   return 0;
 };

src/sphSmooth.hpp

@@ -82,10 +82,14 @@ namespace CosmoTool
       return internal.currentCenter;
     }
 
+    template<typename FuncT>
     ComputePrecision computeSmoothedValue(const typename SPHTree::coords& c,
-					  computeParticleValue fun, SPHState *state = 0);    
+					  FuncT fun, SPHState *state = 0);    
+	
+	template<typename FuncT>
     ComputePrecision computeInterpolatedValue(const typename SPHTree::coords& c,
-					      computeParticleValue fun, SPHState *state = 0);    
+					      FuncT fun, SPHState *state = 0);
+					          
     ComputePrecision getMaxDistance(const typename SPHTree::coords& c, 
 				    SPHNode *node) const;
 
@@ -101,7 +105,8 @@ namespace CosmoTool
     }
     // END
 
-    void runForEachNeighbour(runParticleValue fun, SPHState *state = 0);
+    template<typename FuncT>
+    void runForEachNeighbour(FuncT fun, SPHState *state = 0);
     void addGridSite(const typename SPHTree::coords& c);
 
     bool hasNeighbours() const;
@@ -127,12 +132,15 @@ namespace CosmoTool
     uint32_t maxNgb;
     SPHTree *tree;
     
+    template<typename FuncT>
     ComputePrecision computeWValue(const typename SPHTree::coords & c,
 				   SPHCell& cell,
 				   CoordType d,
-				   computeParticleValue fun, SPHState *state);
+				   FuncT fun, SPHState *state);
+    
+    template<typename FuncT>
     void runUnrollNode(SPHNode *node,
-		       runParticleValue fun);
+		       FuncT fun);
   };
 
   template<class ValType1, class ValType2, int Ndims>

src/sphSmooth.tcc

@@ -21,10 +21,11 @@ SPHSmooth<ValType,Ndims>::~SPHSmooth()
 }
 
 template<typename ValType, int Ndims>
+template<typename FuncT>
 ComputePrecision SPHSmooth<ValType,Ndims>::computeWValue(const typename SPHTree::coords& c,
 							 SPHCell& cell,
 							 CoordType d,
-							 computeParticleValue fun, SPHState *state)
+							 FuncT fun, SPHState *state)
 {
   CoordType weight;
 
@@ -117,9 +118,10 @@ SPHSmooth<ValType,Ndims>::fetchNeighboursOnVolume(const typename SPHTree::coords
 }				
 
 template<typename ValType, int Ndims>
+template<typename FuncT>
 ComputePrecision 
 SPHSmooth<ValType,Ndims>::computeSmoothedValue(const typename SPHTree::coords& c,
-                                               computeParticleValue fun, SPHState *state)
+                                               FuncT fun, SPHState *state)
 {
   if (state == 0)
     state = &internal;
@@ -144,8 +146,9 @@ ComputePrecision interpolateOne(const ValType& t)
 
 // WARNING ! Cell's weight must be 1 !!!
 template<typename ValType, int Ndims>
+template<typename FuncT>
 ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typename SPHTree::coords& c,
-								    computeParticleValue fun, SPHState *state)
+								    FuncT fun, SPHState *state)
 {
   if (state == 0)
     state = &internal;
@@ -164,7 +167,8 @@ ComputePrecision SPHSmooth<ValType,Ndims>::computeInterpolatedValue(const typena
 }
 
 template<typename ValType, int Ndims>
-void SPHSmooth<ValType,Ndims>::runForEachNeighbour(runParticleValue fun, SPHState *state)
+template<typename FuncT>
+void SPHSmooth<ValType,Ndims>::runForEachNeighbour(FuncT fun, SPHState *state)
 {
   if (state == 0)
     state = &internal;