First performance improvements

examples/CMakeLists.txt

@@ -8,3 +8,4 @@ add_subdirectory(point_charge_coulomb)
 add_subdirectory(point_charge_cherenkov)
 add_subdirectory(point_charge_askaryan)
 add_subdirectory(shower_profile)
+add_subdirectory(impacting_shower)

examples/dipole_ice/dipole_ice.cxx

@@ -26,8 +26,6 @@ int main(int argc, char* argv[]) {
   std::string wf_path = argv[1];
 
   auto eps = [](scalar_t r, scalar_t z) {
-
-    return 1.0;
 
     scalar_t z_m = z / 3.0;    
     if(z_m > 0.0) {
@@ -58,14 +56,14 @@ int main(int argc, char* argv[]) {
 
   // CylinderGeometry geom(20, -15, 15, eps);
   // InfEDipoleAntenna dipole(0.0, 10.0, 0.0, impulse_response);
+  // scalar_t t_end = 250;
+
+  CylinderGeometry geom(300, -300, 300, eps);
+  InfEDipoleAntenna dipole(0.0, 10.0, -100.0, impulse_response);
+  scalar_t t_end = 450;
 
-  CylinderGeometry geom(70, -150, 150, eps);
-  InfEDipoleAntenna dipole(0.0, 10.0, 0.0, impulse_response);
-
-  scalar_t t_end = 150;
-  //scalar_t t_end = 25;
   CylindricalWeightingFieldCalculator wfc(geom, dipole, t_end);
-  wfc.Calculate(wf_path);
+  wfc.Calculate(wf_path, "/scratch/midway3/windischhofer/eisvogel/");
 
   return 0;
 }

examples/dipole_ice/run_dipole_ice.sbatch

@@ -0,0 +1,15 @@
+#!/bin/bash
+#SBATCH --job-name=dipole_ice
+#SBATCH --output=/home/windischhofer/data/windischhofer/eisvogel/dipole_ice.out
+#SBATCH --error=/home/windischhofer/data/windischhofer/eisvogel/dipole_ice.err
+#SBATCH --account=pi-avieregg
+#SBATCH --time=24:00:00
+#SBATCH --partition=avieregg
+#SBATCH --nodes=1
+#SBATCH --exclusive
+#SBATCH --ntasks-per-node=64
+#SBATCH --mem=257660
+#SBATCH --mail-type=NONE
+#SBATCH --mail-user=None
+
+sh /home/windischhofer/Eisvogel/examples/dipole_ice/run_dipole_ice.sh

examples/dipole_ice/run_dipole_ice.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+
+export EISVOGELDIR=/home/windischhofer/Eisvogel/
+source ${EISVOGELDIR}/setup_mdwy.sh
+
+cd ${EISVOGELDIR}/build/
+source ./setup.sh
+mpirun --mca orte_base_help_aggregate 0 -np 64 examples/dipole_ice/dipole_ice /home/windischhofer/data/windischhofer/eisvogel/wf_dipole_ice_deep_meep_exp_ratio_20_300

examples/impacting_shower/CMakeLists.txt

@@ -0,0 +1,5 @@
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED YES)
+
+add_executable(impacting_shower impacting_shower.cxx)
+target_link_libraries(impacting_shower eisvogel)

examples/impacting_shower/impacting_shower.cxx

@@ -0,0 +1,65 @@
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cmath>
+
+#include "Eisvogel/Common.hh"
+#include "Eisvogel/SignalCalculator.hh"
+#include "Eisvogel/Current0D.hh"
+#include "Eisvogel/SignalExport.hh"
+
+int main(int argc, char* argv[]) {
+
+  if(argc < 2) {
+    throw;
+  }
+
+  std::string wf_path = argv[1];
+  SignalCalculator calc(wf_path);
+
+  // test trajectory: a point charge moving parallel to the z-axis 
+  // with a constant impact parameter of 'b' along the z-axis
+  scalar_t b = 100;
+  scalar_t zmax = 0.0;
+  scalar_t zscale = 10;
+  scalar_t tstart = -40, tend = 40;
+  scalar_t beta = 1.0;  // for in-ice shower
+  scalar_t qmax = 1e3;
+  auto charge = [&](scalar_t z){
+    scalar_t cur_charge = qmax * std::exp(-std::pow((z - zmax) / zscale, 2.0));
+    if(cur_charge < 1) {
+      cur_charge = 0;
+    }
+    return cur_charge;
+  };
+
+  std::cout << "Building trajectory ..." << std::endl;
+
+  std::vector<CoordVector> points;
+  std::vector<scalar_t> charges;
+
+  for(scalar_t cur_t = tstart; cur_t < tend; cur_t += 0.1) {
+    scalar_t cur_z = -beta * cur_t;
+    scalar_t cur_charge = -charge(cur_z);
+
+    std::cout << "z = " << cur_z << ", q = " << cur_charge << std::endl;
+
+    points.push_back(CoordUtils::MakeCoordVectorTXYZ(cur_t, b, 0, cur_z));
+    charges.push_back(cur_charge);
+  }
+  points.push_back(CoordUtils::MakeCoordVectorTXYZ(tend, b, 0, beta * tend));
+  Current0D track(std::move(points), std::move(charges));
+
+  std::cout << "Computing signal ..." << std::endl;
+  std::vector<scalar_t> signal_times, signal_values;
+  for(scalar_t cur_t = 170; cur_t < 280; cur_t += 0.1) {
+    std::cout << "calculating signal for t = " << cur_t << std::endl;
+    scalar_t cur_signal = calc.ComputeSignal(track, cur_t);
+    signal_times.push_back(cur_t);
+    signal_values.push_back(cur_signal);
+  }
+
+  ExportSignal(signal_times, signal_values, "./impacting_shower_signal.csv");
+
+  return 0;
+}

examples/point_charge_cherenkov/point_charge_cherenkov.cxx

@@ -18,22 +18,23 @@ int main(int argc, char* argv[]) {
   // 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 tstart = -120, tend = -1;
   scalar_t charge = 1;
-  scalar_t beta = 1.2;
+  scalar_t beta = 0.99;
 
   std::cout << "Building trajectory ..." << std::endl;
   Current0D track({
-      CoordUtils::MakeCoordVectorTXYZ(tstart, b, 0, beta * tstart),
-      CoordUtils::MakeCoordVectorTXYZ(tend, b, 0, beta * tend)
+      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 = 10; cur_t < 45; cur_t += 0.1) {
+  //for(scalar_t cur_t = 10; cur_t < 45; cur_t += 0.1) {
+  for(scalar_t cur_t = 40; cur_t < 110; 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);

include/Eisvogel/CylindricalWeightingFieldCalculator.hh

@@ -13,7 +13,7 @@ class CylindricalWeightingFieldCalculator {
 public:
   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 = "");
+  void Calculate(std::filesystem::path outdir, std::filesystem::path mergedir = "", std::filesystem::path tmpdir = "");
 
 private:
 

include/Eisvogel/WeightingField.hh

@@ -60,6 +60,22 @@ public:
 
   template <typename KernelT = DefaultKernel>
   scalar_t E_phi(CoordVector pos);
+
+  // // Vectorized accessors for field components
+  // template <typename KernelT = DefaultKernel>
+  // std::vector<scalar_t> E_r(std::vector<CoordVector>& pos);
+
+  // template <typename KernelT = DefaultKernel>
+  // std::vector<scalar_t> E_z(std::vector<CoordVector>& pos);
+
+  // template <typename KernelT = DefaultKernel>
+  // std::vector<scalar_t> E_phi(std::vector<CoordVector>& pos);
+
+  // template <typename KernelT = DefaultKernel>
+  // scalar_t Eip(CoordVector& pos, CoordVector& vec);
+
+  // template <typename KernelT = DefaultKernel>
+  // std::vector<scalar_t> Eip(std::vector<CoordVector>& pos, std::vector<CoordVector>& vec);
 
   DeltaVector GetSamplingIntervals() const;
   CoordVector GetStartCoords() const;

src/CylindricalWeightingFieldCalculator.cxx

@@ -325,32 +325,42 @@ namespace meep {
 
 } // end namespace meep
 
-void CylindricalWeightingFieldCalculator::Calculate(std::filesystem::path outdir, std::filesystem::path tmpdir) {
+void CylindricalWeightingFieldCalculator::Calculate(std::filesystem::path outdir, std::filesystem::path mergedir, std::filesystem::path tmpdir) {
 
+  // Prepare merge directory
+  if(mergedir.empty()) {
+    mergedir = outdir;
+  }
+
   // TODO: this will get a lot easier once we can have all three components together in the same array
-  std::filesystem::path outdir_Er = outdir / "E_r";
-  std::filesystem::path outdir_Ez = outdir / "E_z";
+  std::filesystem::path mergedir_E_r = mergedir / "E_r";
+  std::filesystem::path mergedir_E_z = mergedir / "E_z";
 
   if(meep::am_master()) {
 
-    // Prepare output directory
-    if(!std::filesystem::exists(outdir)) {
-      std::filesystem::create_directory(outdir);
+    if(!std::filesystem::exists(mergedir)) {
+      std::filesystem::create_directory(mergedir);
     }
 
-    std::filesystem::create_directory(outdir_Er);
-    std::filesystem::create_directory(outdir_Ez);
-  }
+    if(!std::filesystem::exists(mergedir_E_r)) {
+      std::filesystem::create_directory(mergedir_E_r);
+    }
 
-  meep::all_wait();
+    if(!std::filesystem::exists(mergedir_E_z)) {
+      std::filesystem::create_directory(mergedir_E_z);
+    }
+  }
 
   // Make sure to have a valid temporary directory to store the results as we go along
   if(tmpdir.empty()) {
     char tmpdir_template[] = "/tmp/eisvogel.XXXXXX";
-    tmpdir = std::string(mkdtemp(tmpdir_template));
+    tmpdir = std::filesystem::path(mkdtemp(tmpdir_template));
   }
-
+  
   std::cout << "Using tmpdir = " << tmpdir << std::endl;
+  std::cout << "Using mergedir = " << mergedir << std::endl;
+
+  meep::all_wait();
 
   // This is only for cross-checking the geometry for now
   // f -> output_hdf5(meep::Dielectric, gv -> surroundings());
@@ -386,7 +396,11 @@ void CylindricalWeightingFieldCalculator::Calculate(std::filesystem::path outdir
     cld.ind_t = stepcnt++;
     m_f -> loop_in_chunks(meep::eisvogel_saving_chunkloop, static_cast<void*>(&cld), m_f -> total_volume());
 
+    std::cout << "LLL now on stepcnt = " << std::endl;
+    std::cout << stepcnt << std::endl;
+
     if((stepcnt % 400) == 0) {
+      std::cout << "BBB now merging chunks" << std::endl;
       fstor -> MergeChunks(0, 400);
     }
   }
@@ -396,9 +410,9 @@ void CylindricalWeightingFieldCalculator::Calculate(std::filesystem::path outdir
   // TODO: again, will get better once the three separate arrays are gone
   // TODO: for large weighting fields, will have to move chunks to the permanent location continuously throughout the calculation so as not to fill up local storage
   //       move them so that only the complete chunks (surviving after defragmentation) are moved that won't need to be accessed anymore
-  std::cout << "moving output into final location ...";
-  std::filesystem::copy(tmpdir / "E_r", outdir_Er, std::filesystem::copy_options::recursive);
-  std::filesystem::copy(tmpdir / "E_z", outdir_Ez, std::filesystem::copy_options::recursive);
+  std::cout << "moving output into merging location ...";
+  std::filesystem::copy(tmpdir / "E_r", mergedir_E_r, std::filesystem::copy_options::recursive);
+  std::filesystem::copy(tmpdir / "E_z", mergedir_E_z, std::filesystem::copy_options::recursive);
   std::cout << " done!" << std::endl;
 
   // Wait until everybody has finished copying
@@ -411,8 +425,10 @@ void CylindricalWeightingFieldCalculator::Calculate(std::filesystem::path outdir
   std::cout << "==============================================" << std::endl;
 
   if(meep::am_master()) {
-    std::shared_ptr<CylindricalWeightingField> cwf = std::make_shared<CylindricalWeightingField>(outdir, *m_start_coords, *m_end_coords);
+    std::shared_ptr<CylindricalWeightingField> cwf = std::make_shared<CylindricalWeightingField>(mergedir, *m_start_coords, *m_end_coords);
     cwf -> MakeMetadataPersistent();    
     cwf -> RebuildChunks(requested_chunk_size);
+
+    std::filesystem::copy(mergedir, outdir, std::filesystem::copy_options::recursive);
   }
 }

src/DenseNDArray.hh

@@ -157,6 +157,14 @@ public:
     return rhs.m_data == m_data;
   }
 
+  // to assign blocks of data in an efficient way
+  void copy_from(const DenseNDArray<T, dims>& src,
+		 const DenseNDArray<std::size_t, 1>& ind_src_start, const DenseNDArray<std::size_t, 1>& ind_src_stop,
+		 const DenseNDArray<std::size_t, 1>& ind_dest_start, const DenseNDArray<std::size_t, 1>& ind_dest_stop) {
+
+
+  }
+
   // template <std::size_t len>
   // bool operator==(const std::array<T, len>& rhs) const requires(dims == 1) {
   //   if(len != m_data.size()) {