can also store now

include/Eisvogel/SignalCalculator.hh

@@ -17,10 +17,6 @@ public:
   SignalCalculator(const std::string& geometry_path);
   scalar_t ComputeSignal(const Current0D& track, scalar_t t_sig);
   scalar_t ComputeSignal(const SparseCurrentDensity3D& current_distribution, scalar_t t_sig);
-
-public:
-
-  using CylindricalWeightingField = WeightingField<TRZFieldIndexer<CylindricalSymmetry>, RZFieldStorage>;
 
 private:
   std::string m_geometry_path;

include/Eisvogel/WeightingField.hh

@@ -18,13 +18,19 @@ class TRZFieldIndexer {
 public:
 
   TRZFieldIndexer(std::filesystem::path index_path, FieldStorage& storage);
+  TRZFieldIndexer(std::filesystem::path index_path, const CoordVector& start_coords, const CoordVector& end_coords);
 
   CoordVector GetFieldIndexFromCoord(CoordVector pos_txyz);
   IndexVector GetFieldIndex(GridVector inds_trz);
 
   DeltaVector GetSamplingIntervals() const;
 
+  void MakePersistent();
+
 private:
+
+  std::filesystem::path m_meta_path;
+
   std::shared_ptr<IndexVector> m_shape;
   std::shared_ptr<CoordVector> m_start_coords;
   std::shared_ptr<CoordVector> m_end_coords;  
@@ -36,7 +42,11 @@ class WeightingField {
 public:
 
   WeightingField(std::string wf_path);
+  WeightingField(std::string wf_path, const CoordVector& start_coords, const CoordVector& end_coords);
 
+  template <typename ...Params>
+  void RegisterChunk(Params&&... params);
+
   // Accessors for field components
   template <typename KernelT>
   scalar_t E_r(CoordVector pos);
@@ -48,6 +58,8 @@ public:
   scalar_t E_phi(CoordVector pos);
 
   DeltaVector GetSamplingIntervals() const;
+
+  void MakeMetadataPersistent();
 
 private:
 
@@ -65,4 +77,7 @@ private:
 
 #include "WeightingField.hxx"
 
+#include "Symmetry.hh"
+using CylindricalWeightingField = WeightingField<TRZFieldIndexer<CylindricalSymmetry>, RZFieldStorage>;
+
 #endif

src/FieldStorage.cxx

@@ -25,6 +25,16 @@ scalar_t RZFieldStorage::E_phi(IndexVector& ind) {
   return 0.0;
 }
 
+void RZFieldStorage::RegisterChunk(const chunk_t& chunk_E_r, const chunk_t& chunk_E_z, const IndexVector chunk_start_inds) {
+  m_E_r -> RegisterChunk(chunk_E_r, chunk_start_inds);
+  m_E_z -> RegisterChunk(chunk_E_z, chunk_start_inds);
+}
+
+void RZFieldStorage::MakeIndexPersistent() {
+  m_E_r -> FlushIndex();
+  m_E_z -> FlushIndex();
+}
+
 // calculate from the stored components
 // scalar E_x(IndexVector& ind) { ... };
 // scalar E_y(IndexVector& ind) { ... };

src/FieldStorage.hh

@@ -15,6 +15,7 @@ class RZFieldStorage : public FieldStorage {
 
 public:
   using storage_t = DistributedNDArray<scalar_t, 3>;
+  using chunk_t = DenseNDArray<scalar_t, 3>;
 
 public:
 
@@ -26,6 +27,10 @@ public:
 
   IndexVector shape();
 
+  void RegisterChunk(const chunk_t& chunk_E_r, const chunk_t& chunk_E_z, const IndexVector chunk_start_inds);
+
+  void MakeIndexPersistent();
+
   // vector E_rzphi(IndexVector& ind);
 
   // calculate from the stored components

src/SignalCalculator.cxx

@@ -22,10 +22,10 @@ scalar_t SignalCalculator::ComputeSignal(const Current0D& track, scalar_t t_sig)
   // the shower (still need to pay attention that we're not generating I/O unecessarily)
 
   // Compute the signal
-  scalar_t sig_old = m_integrator.integrate(t_sig, track);
+  // scalar_t sig_old = m_integrator.integrate(t_sig, track);
   scalar_t sig_new = integrate<CylindricalWeightingField>(*m_wf, t_sig, track);
 
-  std::cout << "old = " << sig_old << "  < --- >   new = " << sig_new << std::endl;
+  // std::cout << "old = " << sig_old << "  < --- >   new = " << sig_new << std::endl;
 
   return sig_new;
 }

src/WeightingField.hxx

@@ -1,17 +1,35 @@
 #include "Eisvogel/Serialization.hh"
 
 template <typename SymmetryT>
-TRZFieldIndexer<SymmetryT>::TRZFieldIndexer(std::filesystem::path index_path, FieldStorage& storage) {
+TRZFieldIndexer<SymmetryT>::TRZFieldIndexer(std::filesystem::path index_path, FieldStorage& storage) :
+  m_meta_path(index_path / "wf_meta.bin") {
 
   std::fstream ifs;
-  std::filesystem::path meta_path = index_path / "wf_meta.bin";
-  ifs.open(meta_path, std::ios::in | std::ios::binary);
+  ifs.open(m_meta_path, std::ios::in | std::ios::binary);
   stor::Serializer iser(ifs);
   m_start_coords = std::make_shared<CoordVector>(iser.deserialize<CoordVector>());
   m_end_coords = std::make_shared<CoordVector>(iser.deserialize<CoordVector>());
   m_shape = std::make_shared<IndexVector>(storage.shape());
 }
 
+template <typename SymmetryT>
+TRZFieldIndexer<SymmetryT>::TRZFieldIndexer(std::filesystem::path index_path, const CoordVector& start_coords, const CoordVector& end_coords) :
+  m_meta_path(index_path / "wf_meta.bin") {
+
+  m_start_coords = std::make_shared<CoordVector>(start_coords);
+  m_end_coords = std::make_shared<CoordVector>(end_coords);
+}
+
+template <typename SymmetryT>
+void TRZFieldIndexer<SymmetryT>::MakePersistent() {
+
+  std::fstream ofs;
+  ofs.open(m_meta_path, std::ios::out | std::ios::binary);  
+  stor::Serializer oser(ofs);
+  oser.serialize(*m_start_coords);
+  oser.serialize(*m_end_coords);
+}
+
 template <typename SymmetryT>
 CoordVector TRZFieldIndexer<SymmetryT>::GetSamplingIntervals() const {
   return (*m_end_coords - *m_start_coords) / (CoordVector)*m_shape;
@@ -42,13 +60,43 @@ IndexVector TRZFieldIndexer<SymmetryT>::GetFieldIndex(GridVector inds_trz) {
 
 // ---------------
 
+template <class FieldIndexerT, class FieldStorageT>
+WeightingField<FieldIndexerT, FieldStorageT>::WeightingField(std::string wf_path, const CoordVector& start_coords, const CoordVector& end_coords) :
+  m_wf_path(wf_path) {
+
+  // we might have to create the directory structure if this is an empty weighting field
+  if(!std::filesystem::exists(m_wf_path)) {
+    std::filesystem::create_directory(m_wf_path);
+  }
+
+  m_field_storage = std::make_shared<FieldStorageT>(wf_path, 10);
+  m_field_indexer = std::make_shared<FieldIndexerT>(wf_path, start_coords, end_coords);  
+}
+
 template <class FieldIndexerT, class FieldStorageT>
 WeightingField<FieldIndexerT, FieldStorageT>::WeightingField(std::string wf_path) : m_wf_path(wf_path) {
 
+  // we might have to create the directory structure if this is an empty weighting field
+  if(!std::filesystem::exists(m_wf_path)) {
+    throw std::runtime_error("No usable weighting field found at this location!");
+  }
+
   m_field_storage = std::make_shared<FieldStorageT>(wf_path, 10);
   m_field_indexer = std::make_shared<FieldIndexerT>(wf_path, *m_field_storage);  
 }
 
+template <class FieldIndexerT, class FieldStorageT>
+template <typename ...Params>
+void WeightingField<FieldIndexerT, FieldStorageT>::RegisterChunk(Params&&... params) {
+  m_field_storage -> RegisterChunk(std::forward<Params>(params)...);
+}
+
+template <class FieldIndexerT, class FieldStorageT>
+void WeightingField<FieldIndexerT, FieldStorageT>::MakeMetadataPersistent() {
+  m_field_storage -> MakeIndexPersistent();
+  m_field_indexer -> MakePersistent();
+}
+
 template <class FieldIndexerT, class FieldStorageT>
 DeltaVector WeightingField<FieldIndexerT, FieldStorageT>::GetSamplingIntervals() const {
   return m_field_indexer -> GetSamplingIntervals();

src/WeightingFieldUtils.cxx

@@ -5,6 +5,7 @@
 #include "Eisvogel/MathUtils.hh"
 #include "Eisvogel/Serialization.hh"
 #include "Eisvogel/DistributedWeightingField.hh"
+#include "Eisvogel/WeightingField.hh"
 #include <cmath>
 #include <iostream>
 #include <filesystem>
@@ -21,7 +22,8 @@ namespace WeightingFieldUtils {
     // make sure to start from scratch
     std::filesystem::remove_all(wf_path);
 
-    DistributedWeightingField dwf(wf_path, start_coords, end_coords);
+    //DistributedWeightingField dwf(wf_path, start_coords, end_coords);
+    CylindricalWeightingField cwf(wf_path, start_coords, end_coords);
 
     // compute required step size for sampling of weighting field
     scalar_t c = 1.0;  // speed of light in vacuum
@@ -67,11 +69,13 @@ namespace WeightingFieldUtils {
 					      tp, N, r_min, os_factor, n);
 
       // Register chunk buffers
-      dwf.RegisterChunk(chunk_buffer_E_r, chunk_buffer_E_z, chunk_buffer_E_phi, chunk_start_inds);
+      // dwf.RegisterChunk(chunk_buffer_E_r, chunk_buffer_E_z, chunk_buffer_E_phi, chunk_start_inds);
+      cwf.RegisterChunk(chunk_buffer_E_r, chunk_buffer_E_z, chunk_start_inds);
 
     }   
 
-    dwf.Flush();
+    // dwf.Flush();
+    cwf.MakeMetadataPersistent();
   }
 
   // TODO: three separate `ScalarField3D`s to be replaced with single vector field