Sparsification and storage (#25)

* Handles storage of dense and sparse field chunks

* Includes functionality to merge chunks and re-chunk the full field for faster memory access

examples/dipole_free_space/dipole_free_space.cxx

@@ -13,15 +13,15 @@ int main(int argc, char* argv[]) {
   std::string wf_path = argv[1];
 
   // Domain of weighting field
-  CoordVector start_coords = CU::MakeCoordVectorTRZ(0.0, 0.0, -20.0);
-  CoordVector end_coords = CU::MakeCoordVectorTRZ(25.0, 20.0, 20.0);
-
+  CoordVector start_coords = CU::MakeCoordVectorTRZ(0.0, 0.0, -150.0);
+  CoordVector end_coords = CU::MakeCoordVectorTRZ(150.0, 70.0, 150.0);
+  
   // Filter parameters
   scalar_t tp = 2.0;
   unsigned int N = 4;
 
   // Sampling parameters
-  scalar_t os_factor = 50;
+  scalar_t os_factor = 30;
   scalar_t r_min = 0.1;
 
   WFU::CreateElectricDipoleWeightingField(wf_path, start_coords, end_coords, tp, N, r_min, os_factor);

examples/dipole_ice/dipole_ice.cxx

@@ -1,7 +1,7 @@
 #include <iostream>
 #include <cmath>
 #include "Eisvogel/Common.hh"
-#include "Eisvogel/WeightingFieldCalculator.hh"
+#include "Eisvogel/CylindricalWeightingFieldCalculator.hh"
 #include "Eisvogel/Antenna.hh"
 #include "Eisvogel/CoordUtils.hh"
 
@@ -26,7 +26,25 @@ int main(int argc, char* argv[]) {
   std::string wf_path = argv[1];
 
   auto eps = [](scalar_t r, scalar_t z) {
-    return 1.00;
+
+    return 1.0;
+
+    scalar_t z_m = z / 3.0;    
+    if(z_m > 0.0) {
+      return 1.0;
+    }
+
+    scalar_t density = 0.0;
+
+    if(z_m > -14.9) {
+      density = 0.917 - 0.594 * std::exp(z_m / 30.8);
+    }
+    else {
+      density = 0.917 - 0.367 * std::exp((z_m + 14.9) / 40.5);
+    }
+
+    double eps = std::pow(1 + 0.845 * density, 2.0);
+    return eps;
   };
 
   auto impulse_response = [](scalar_t t) {
@@ -37,12 +55,16 @@ int main(int argc, char* argv[]) {
     }
     return std::pow(t / tp * N, N) * std::exp(-t / tp * N) / (tp * std::exp(std::lgamma(N)));
   };
+
+  // CylinderGeometry geom(20, -15, 15, eps);
+  // InfEDipoleAntenna dipole(0.0, 10.0, 0.0, impulse_response);
 
-  CylinderGeometry geom(20, -20, 20, eps);
-  InfEDipoleAntenna dipole(0.0, 10.0, -2.0, impulse_response);
+  CylinderGeometry geom(70, -150, 150, eps);
+  InfEDipoleAntenna dipole(0.0, 10.0, 0.0, impulse_response);
 
-  scalar_t t_end = 25;
-  WeightingFieldCalculator wfc(geom, dipole, t_end);
+  scalar_t t_end = 150;
+  //scalar_t t_end = 25;
+  CylindricalWeightingFieldCalculator wfc(geom, dipole, t_end);
   wfc.Calculate(wf_path);
 
   return 0;

examples/plot_weighting_field/plot_weighting_field.py

@@ -20,7 +20,7 @@ def plot_2d(outpath, vals_xy, vals_z, xlabel = "", ylabel = "", zlabel = "", tit
     vals_mesh_x, vals_mesh_y = np.meshgrid(vals_fine_x, vals_fine_y)
     vals_interp_z = interpolator(vals_mesh_x, vals_mesh_y)
 
-    fig, ax = plt.subplots(1, 1, figsize = (5, 5))
+    fig, ax = plt.subplots(1, 1)
     surveyplot = ax.pcolormesh(vals_mesh_x, vals_mesh_y, vals_interp_z, cmap = cmap,
                                norm = matplotlib.colors.SymLogNorm(linthresh = 1e-6))
     divider = make_axes_locatable(ax)
@@ -36,6 +36,8 @@ def plot_2d(outpath, vals_xy, vals_z, xlabel = "", ylabel = "", zlabel = "", tit
     ax.set_ylabel(ylabel, fontsize = fs)
     ax.set_title(title, fontsize = fs)
 
+    ax.set_aspect('equal')
+
     plt.tight_layout()
 
     fig.savefig(outpath, dpi = 300)
@@ -50,8 +52,8 @@ def plot_fields(outdir, cwf, tval, num_pts = 200):
     vals_E_r = []
     vals_E_z = []
     vals_E_abs = []
-    for cur_xval in np.linspace(start_coords_txyz[1] + 2, end_coords_txyz[1] - 2, num_pts):
-        for cur_zval in np.linspace(start_coords_txyz[3] + 2, end_coords_txyz[3] - 2, num_pts):
+    for cur_xval in np.linspace(start_coords_txyz[1] + 2, end_coords_txyz[1] - 2, int(end_coords_txyz[1] - start_coords_txyz[1])):
+        for cur_zval in np.linspace(start_coords_txyz[3] + 2, end_coords_txyz[3] - 2, int(end_coords_txyz[3] - start_coords_txyz[3])):
 
             vals_xz.append([cur_xval, cur_zval])
             cur_yval = 0.0

examples/point_charge_cherenkov/point_charge_cherenkov.cxx

@@ -15,25 +15,25 @@ int main(int argc, char* argv[]) {
   std::string wf_path = argv[1];
   SignalCalculator calc(wf_path);
 
-  // test trajectory: a point charge moving parallel to the x-axis 
-  // with a constant impact parameter of 'b' along the z-axis
-  scalar_t b = 20;
-  scalar_t tstart = -250, tend = 250;
+  // test trajectory: a point charge moving parallel to the z-axis 
+  // with a constant impact parameter of 'b' along the x-axis
+  scalar_t b = 50;
+  scalar_t tstart = -100, tend = 100;
   scalar_t charge = 1;
   scalar_t beta = 1.2;
 
   std::cout << "Building trajectory ..." << std::endl;
   Current0D track({
-      CoordUtils::MakeCoordVectorTXYZ(tstart, beta * tstart, 0, b),
-  	CoordUtils::MakeCoordVectorTXYZ(tend, beta * tend, 0, b)
-  	},
+      CoordUtils::MakeCoordVectorTXYZ(tstart, b, 0, beta * tstart),
+      CoordUtils::MakeCoordVectorTXYZ(tend, b, 0, beta * tend)
+    },
     {charge}
     );
 
   std::cout << "Computing signal ..." << std::endl;
 
   std::vector<scalar_t> signal_times, signal_values;
-  for(scalar_t cur_t = 0; cur_t < 40; cur_t += 1) {
+  for(scalar_t cur_t = 10; cur_t < 45; cur_t += 0.1) {
     scalar_t cur_signal = calc.ComputeSignal(track, cur_t);
     signal_times.push_back(cur_t);
     signal_values.push_back(cur_signal);

extern/shower_profile/CMakeLists.txt

@@ -18,12 +18,14 @@ target_include_directories(Shower1D PUBLIC
                           "${PROJECT_SOURCE_DIR}/extern/utilities"
                           "${PROJECT_SOURCE_DIR}/extern/constants"
                           "${PROJECT_SOURCE_DIR}/include/"
+			  "${PROJECT_SOURCE_DIR}/src/"
                           )
 
 target_include_directories(ShowerCreator PUBLIC
                           "${PROJECT_SOURCE_DIR}/extern/shower_profile/shower_1D"
                           "${PROJECT_SOURCE_DIR}/extern/shower_profile/shower_2D"
                           "${PROJECT_SOURCE_DIR}/extern/utilities"
+			  "${PROJECT_SOURCE_DIR}/src/"
                           )
 target_include_directories(Environment PUBLIC
                           "${PROJECT_SOURCE_DIR}/extern/utilities"
@@ -35,4 +37,5 @@ target_include_directories(Shower2D PUBLIC
                         "${PROJECT_SOURCE_DIR}/extern/utilities"
                         "${PROJECT_SOURCE_DIR}/extern/constants"
                         "${PROJECT_SOURCE_DIR}/include/"
-                        )
+			"${PROJECT_SOURCE_DIR}/src/"
+                        )

extern/shower_profile/charge_excess_profile/charge_excess_profile.cpp

@@ -2,7 +2,7 @@
 #define CHARGEEXCESSPROFILE_CLASS
 #include <vector>
 #include "Eisvogel/Common.hh"
-#include "Eisvogel/NDArray.hh"
+#include "DenseNDArray.hh"
 
 namespace showers {
 class ChargeExcessProfile {

extern/shower_profile/shower_creator/shower_creator.h

@@ -4,7 +4,7 @@
 #include "shower_2D.h"
 #include "ice_profile.h"
 #include "charge_excess_profile.cpp"
-#include "Eisvogel/NDArray.hh"
+#include "NDArray.hh"
 #include <array>
 #include <string>
 #include <map>

include/Eisvogel/CoordUtils.hh

@@ -2,7 +2,7 @@
 #define __COORD_UTILS_HH
 
 #include "Common.hh"
-#include "NDArray.hh"
+#include "DenseNDArray.hh"
 #include <cmath>
 
 using CoordVector = DenseVector<scalar_t>;

include/Eisvogel/WeightingFieldCalculator.hh → include/Eisvogel/CylindricalWeightingFieldCalculator.hh

@@ -1,5 +1,5 @@
-#ifndef __WEIGHTING_FIELD_CALCULATOR__HH
-#define __WEIGHTING_FIELD_CALCULATOR__HH
+#ifndef __CYLINDRICAL_WEIGHTING_FIELD_CALCULATOR__HH
+#define __CYLINDRICAL_WEIGHTING_FIELD_CALCULATOR__HH
 
 #include <memory>
 #include <string>
@@ -8,11 +8,11 @@
 #include "Antenna.hh"
 #include "Geometry.hh"
 
-class WeightingFieldCalculator {
+class CylindricalWeightingFieldCalculator {
 
 public:
-  WeightingFieldCalculator(CylinderGeometry& geom, const Antenna& antenna, scalar_t t_end,
-			   double courant_factor = 0.5, double resolution = 40, double pml_width = 1.0);
+  CylindricalWeightingFieldCalculator(CylinderGeometry& geom, const Antenna& antenna, scalar_t t_end,
+				      double courant_factor = 0.5, double resolution = 12, double pml_width = 1.0);
   void Calculate(std::filesystem::path outdir, std::filesystem::path tmpdir = "");
 
 private:

include/Eisvogel/Geometry.hh

@@ -31,6 +31,10 @@ public:
   meep::vec toMeepCoords(const CoordVector& coords) {
     return meep::veccyl(CoordUtils::getR(coords), CoordUtils::getZ(coords) - z_min);
   };
+
+  CoordVector toEisvogelCoords(const meep::vec& meep_coords) {
+    return CoordUtils::MakeCoordVectorTRZ(0.0, meep_coords.r(), meep_coords.z() + z_min);
+  };
 
 public:
   virtual double eps(const meep::vec& pos);

include/Eisvogel/IteratorUtils.hh

@@ -1,7 +1,8 @@
 #ifndef __ITERATOR_UTILS_HH
 #define __ITERATOR_UTILS_HH
 
-#include "NDArray.hh"
+#include <cuchar>
+#include "DenseNDArray.hh"
 
 // Iterator for dynamic (i.e. known at run time) number of dimensions
 template <typename VecT> class VectorCounter {
@@ -16,7 +17,7 @@ public:
   VectorCounter(const VecT& start, const VecT& end) : m_start(start), m_end(end), m_cur(start), 
 						      number_dims(start.size()) { }
   VectorCounter& operator++() {
-    for(size_t i = 0; i < number_dims; i++) {
+    for(std::size_t i = 0; i < number_dims; i++) {
       if(++m_cur(i) == m_end(i)) {
 	if(i < number_dims - 1) {
 	  m_cur(i) = m_start(i);
@@ -33,7 +34,7 @@ public:
     return m_cur(number_dims - 1) < m_end(number_dims - 1);
   }
 
-  VecT::type& operator()(size_t ind) {return m_cur(ind);}
+  VecT::type& operator()(std::size_t ind) {return m_cur(ind);}
 
   auto begin() {return m_cur.begin();}
   auto end() {return m_cur.end();}

include/Eisvogel/WeightingField.hh

@@ -66,6 +66,9 @@ public:
   scalar_t GetStartCoords(std::size_t dim) const;
   CoordVector GetEndCoords() const;
   scalar_t GetEndCoords(std::size_t dim) const;
+
+  void RebuildChunks(const IndexVector& requested_chunk_size);
+  void MergeChunks(std::size_t dim_to_merge, std::size_t max_dimsize);
 
 private:
 

src/CMakeLists.txt

@@ -25,7 +25,7 @@ if(BUILD_WFC)
     pkg_check_modules(MEEP REQUIRED meep)
 
     add_library(eisvogel_wfc SHARED
-                WeightingFieldCalculator.cxx
+                CylindricalWeightingFieldCalculator.cxx
                 Geometry.cxx
                 Antenna.cxx)
     target_link_directories(eisvogel_wfc PUBLIC ${MEEP_LIBRARY_DIRS})