Development

.gitmodules

@@ -1,3 +1,3 @@
-[submodule "extra/rapidcheck"]
-	path = extras/rapidcheck
+[submodule "rapidcheck"]
+	path = rapidcheck
 	url = https://github.com/emil-e/rapidcheck.git

.travis.yml

@@ -0,0 +1,24 @@
+language: cpp
+compiler:
+  #- clang
+  - gcc
+
+before_script:
+  - sudo ./install/get_boost.sh > /dev/null
+  - sudo ./install/get_eigen.sh > /dev/null
+  - mkdir build && cd build
+  - cmake ..
+
+script: make qtest-bin && ./qtest
+
+addons:
+  apt:
+    sources:
+      - george-edison55-precise-backports
+    packages:
+      - cmake-data
+      - cmake
+      - libfftw3-dev
+        #- libeigen3-dev
+        #- libboost-test-dev
+        #- libboost-program-options-dev

README.md

@@ -1,4 +1,5 @@
-[![DOI](https://zenodo.org/badge/37320094.svg)](https://zenodo.org/badge/latestdoi/37320094)
+[![Build Status](https://travis-ci.org/tbertus/QuEST.svg?branch=development)](https://travis-ci.org/tbertus/QuEST/branch/development) [![DOI](https://zenodo.org/badge/37320094.svg)](https://zenodo.org/badge/latestdoi/37320094)
+
 [![Article](https://img.shields.io/badge/article-Phys.%20Rev.%20A%2096%2C%20033816-B10DC9.svg)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.96.033816)
 
 # QuEST: The Quantum Electromagnetics Simulation Toolbox

install/get_boost.sh

@@ -0,0 +1,14 @@
+#!/bin/bash
+echo "Changing to tmp"
+cd /tmp
+
+echo "Downloading and unpacking boost"
+wget https://dl.bintray.com/boostorg/release/1.68.0/source/boost_1_68_0.tar.bz2 -O boost.tar.bz2
+mkdir -p boost && tar -xvf boost.tar.bz2 -C boost --strip-components 1 > /dev/null
+ls
+
+echo "Building boost"
+cd boost/
+ls
+./bootstrap.sh --with-libraries=program_options,test
+./b2 install

install/get_eigen.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+cd /tmp
+wget https://bitbucket.org/eigen/eigen/get/3.3.2.tar.gz -O eigen.tar.gz
+mkdir -p eigen/build && tar -xvzf eigen.tar.gz -C eigen --strip-components 1
+cd eigen/build && cmake .. && make install

src/integrator/RHS/bloch_rhs.cpp

@@ -12,13 +12,13 @@ Integrator::BlochRHS::BlochRHS(
 {
 }
 
-void Integrator::BlochRHS::evaluate(const int step) const
+void Integrator::BlochRHS::evaluate_present(const int step) const
 {
-  auto eval_and_sum =
-      [step](const InteractionBase::ResultArray &r,
-             const std::shared_ptr<InteractionBase> &interaction) {
-        return r + interaction->evaluate(step);
-      };
+  auto eval_and_sum = [step](
+                          const InteractionBase::ResultArray &r,
+                          const std::shared_ptr<InteractionBase> &interaction) {
+    return r + interaction->evaluate_present_field(step);
+  };
   auto nil = InteractionBase::ResultArray::Zero(num_solutions, 1).eval();
 
   auto projected_efields = std::accumulate(
@@ -28,5 +28,22 @@ void Integrator::BlochRHS::evaluate(const int step) const
     history->array_[solution][step][1] = rhs_functions[solution](
         history->array_[solution][step][0], projected_efields[solution]);
   }
+}
 
+void Integrator::BlochRHS::evaluate(const int step) const
+{
+  auto eval_and_sum = [step](
+                          const InteractionBase::ResultArray &r,
+                          const std::shared_ptr<InteractionBase> &interaction) {
+    return r + interaction->evaluate(step);
+  };
+  auto nil = InteractionBase::ResultArray::Zero(num_solutions, 1).eval();
+
+  auto projected_efields = std::accumulate(
+      interactions.begin(), interactions.end(), nil, eval_and_sum);
+
+  for(int solution = 0; solution < num_solutions; ++solution) {
+    history->array_[solution][step][1] = rhs_functions[solution](
+        history->array_[solution][step][0], projected_efields[solution]);
+  }
 }

src/integrator/RHS/bloch_rhs.h

@@ -17,6 +17,7 @@ class Integrator::BlochRHS : public Integrator::RHS<Eigen::Vector2cd> {
            const std::shared_ptr<History<Eigen::Vector2cd>>,
            std::vector<std::shared_ptr<InteractionBase>>,
            std::vector<BlochFunctionType>);
+  void evaluate_present(const int) const override;
   void evaluate(const int) const override;
 
  private:

src/integrator/RHS/ode_rhs.h

@@ -12,8 +12,10 @@ namespace Integrator {
 
 class Integrator::ODE_RHS : public RHS<double> {
  public:
-  ODE_RHS(const double, const std::shared_ptr<Integrator::History<double>>,
+  ODE_RHS(const double,
+          const std::shared_ptr<Integrator::History<double>>,
           const std::vector<std::function<double(double, double)>>);
+  void evaluate_present(const int) const;
   void evaluate(const int) const;
 
  private:
@@ -36,4 +38,6 @@ void Integrator::ODE_RHS::evaluate(const int n) const
   }
 }
 
+void Integrator::ODE_RHS::evaluate_present(const int n) const { evaluate(n); }
+
 #endif

src/integrator/RHS/rhs.h

@@ -14,6 +14,7 @@ class Integrator::RHS {
  public:
   RHS(const double dt, const std::shared_ptr<History<soltype>> history)
       : dt(dt), history(std::move(history)){};
+  virtual void evaluate_present(const int) const = 0;
   virtual void evaluate(const int) const = 0;
 
  protected:

src/integrator/integrator.h

@@ -80,7 +80,7 @@ void Integrator::PredictorCorrector<soltype>::solve_step(const int step) const
 
   for(int m = 0; m < num_corrector_steps; ++m) {
     corrector(step);
-    rhs->evaluate(step);
+    rhs->evaluate_present(step);
   }
 }
 

src/interactions/AIM/aim_interaction.h

@@ -60,12 +60,20 @@ class AIM::Interaction final : public InteractionBase {
 
   const ResultArray &evaluate(const int t)
   {
-    // I DON'T KNOW WHY THAT NEEDS A CONJUGATE!!!
     results =
         (ff.evaluate(t).conjugate() - nf.evaluate(t)) + direct.evaluate(t);
     return results;
   }
 
+  const ResultArray &evaluate_present_field(const int t)
+  {
+    // I DON'T KNOW WHY THAT NEEDS A CONJUGATE!!!
+    results = (ff.evaluate_present_field(t).conjugate() -
+               nf.evaluate_present_field(t)) +
+              direct.evaluate_present_field(t);
+    return results;
+  }
+
  private:
   std::shared_ptr<Grid> grid;
   std::shared_ptr<Expansions::ExpansionTable> expansion_table;

src/interactions/AIM/direct.cpp

@@ -20,25 +20,58 @@ AIM::DirectInteraction::DirectInteraction(
 const InteractionBase::ResultArray &AIM::DirectInteraction::evaluate(
     const int time_idx)
 {
+  constexpr int RHO_01 = 1;
+
   results.setZero();
+  past_terms_of_convolution.setZero();
 
   for(int pair_idx = 0; pair_idx < shape_[0]; ++pair_idx) {
     const auto &pair = (*interaction_pairs_)[pair_idx];
 
-    for(int i = 0; i < shape_[1]; ++i) {
+    for(int i = 1; i < shape_[1]; ++i) {
       const int s =
           std::max(time_idx - floor_delays_[pair_idx] - i,
                    static_cast<int>(history->array_.index_bases()[1]));
 
-      constexpr int RHO_01 = 1;
-
-      results[pair.first] += (history->array_[pair.second][s][0])[RHO_01] *
-                             coefficients_[pair_idx][i];
-      results[pair.second] += (history->array_[pair.first][s][0])[RHO_01] *
-                              coefficients_[pair_idx][i];
+      past_terms_of_convolution[pair.first] +=
+          (history->array_[pair.second][s][0])[RHO_01] *
+          coefficients_[pair_idx][i];
+      past_terms_of_convolution[pair.second] +=
+          (history->array_[pair.first][s][0])[RHO_01] *
+          coefficients_[pair_idx][i];
     }
+
+    const int s = std::max(time_idx - floor_delays_[pair_idx],
+                           static_cast<int>(history->array_.index_bases()[1]));
+
+    results[pair.first] += (history->array_[pair.second][s][0])[RHO_01] *
+                           coefficients_[pair_idx][0];
+    results[pair.second] += (history->array_[pair.first][s][0])[RHO_01] *
+                            coefficients_[pair_idx][0];
   }
+  results += past_terms_of_convolution;
+  return results;
+}
+
+const InteractionBase::ResultArray &
+AIM::DirectInteraction::evaluate_present_field(const int time_idx)
+{
+  constexpr int RHO_01 = 1;
+
+  results.setZero();
+
+  for(int pair_idx = 0; pair_idx < shape_[0]; ++pair_idx) {
+    const auto &pair = (*interaction_pairs_)[pair_idx];
 
+    const int s = std::max(time_idx - floor_delays_[pair_idx],
+                           static_cast<int>(history->array_.index_bases()[1]));
+
+    results[pair.first] += (history->array_[pair.second][s][0])[RHO_01] *
+                           coefficients_[pair_idx][0];
+    results[pair.second] += (history->array_[pair.first][s][0])[RHO_01] *
+                            coefficients_[pair_idx][0];
+  }
+  results += past_terms_of_convolution;
   return results;
 }
 

src/interactions/AIM/direct.h

@@ -19,6 +19,7 @@ class AIM::DirectInteraction final : public HistoryInteraction {
       std::shared_ptr<const std::vector<Grid::ipair_t>>);
 
   const ResultArray &evaluate(const int) final;
+  const ResultArray &evaluate_present_field(const int) final;
   boost::multi_array<cmplx, 2> coefficient_table(Propagation::Kernel<cmplx> &,
                                                  std::vector<int> &) const;
 

src/interactions/AIM/farfield.cpp

@@ -76,18 +76,59 @@ void AIM::Farfield::propagate(const int step)
 
   Eigen::Map<Eigen::Array3Xcd> observers(&obs_table_(front), 3, nb);
   observers = 0;
+  // initialize temp_observers with the correct size and zero it
+  temp_observers = Eigen::Array3Xcd::Zero(3, nb);
 
-  for(int i = 0; i < table_dimensions_[0]; ++i) {
+  // compute convolution up to but not including the current timestep
+  for(int i = 1; i < table_dimensions_[0]; ++i) {
     // If (step - i) runs "off the end", just propagate src[0][...]
     auto wrap = std::max(step - i, 0) % table_dimensions_[0];
 
     Eigen::Map<Eigen::ArrayXcd> prop(&propagation_table_[i][0][0][0], nb);
     Eigen::Map<Eigen::Array3Xcd> src(&source_table_[wrap][0][0][0][0], 3, nb);
 
     // Use broadcasting to do the x, y, and z component propagation
-    observers += src.rowwise() * prop.transpose();
+    temp_observers += src.rowwise() * prop.transpose();
   }
 
+  // set observers equal to temp_observers
+  // std::cout << temp_observers << "\n\n\n\n" << std::endl;
+  observers += temp_observers;
+
+  auto wrap = std::max(step, 0) % table_dimensions_[0];  // needed ?
+
+  Eigen::Map<Eigen::ArrayXcd> prop(&propagation_table_[0][0][0][0], nb);
+  Eigen::Map<Eigen::Array3Xcd> src(&source_table_[wrap][0][0][0][0], 3, nb);
+  observers += src.rowwise() * prop.transpose();
+
+  const auto o_ptr = &obs_table_(front);
+  fftw_execute_dft(spatial_vector_transforms_.backward,
+                   reinterpret_cast<fftw_complex *>(o_ptr),
+                   reinterpret_cast<fftw_complex *>(o_ptr));
+}
+
+void AIM::Farfield::propagate_present_field(const int step)
+{
+  const auto wrapped_step = step % table_dimensions_[0];
+  const auto nb = 8 * grid->size();
+  const std::array<int, 5> front = {{wrapped_step, 0, 0, 0, 0}};
+
+  const auto s_ptr = &source_table_(front);
+  fftw_execute_dft(spatial_vector_transforms_.forward,
+                   reinterpret_cast<fftw_complex *>(s_ptr),
+                   reinterpret_cast<fftw_complex *>(s_ptr));
+
+  Eigen::Map<Eigen::Array3Xcd> observers(&obs_table_(front), 3, nb);
+  observers = 0;
+
+  observers += temp_observers;
+
+  auto wrap = std::max(step, 0) % table_dimensions_[0];  // needed ?
+
+  Eigen::Map<Eigen::ArrayXcd> prop(&propagation_table_[0][0][0][0], nb);
+  Eigen::Map<Eigen::Array3Xcd> src(&source_table_[wrap][0][0][0][0], 3, nb);
+  observers += src.rowwise() * prop.transpose();
+
   const auto o_ptr = &obs_table_(front);
   fftw_execute_dft(spatial_vector_transforms_.backward,
                    reinterpret_cast<fftw_complex *>(o_ptr),

src/interactions/AIM/farfield.h

@@ -21,12 +21,19 @@ class AIM::Farfield final : public AimBase {
            Normalization::SpatialNorm);
   ~Farfield() = default;
 
+  const ResultArray &evaluate_present_field(const int step) final
+  {
+    fill_source_table(step);
+    propagate_present_field(step);
+    fill_results_table(step);
+
+    return results;
+  }
   const ResultArray &evaluate(const int step) final
   {
     fill_source_table(step);
     propagate(step);
     fill_results_table(step);
-
     return results;
   }
 
@@ -39,11 +46,14 @@ class AIM::Farfield final : public AimBase {
   spacetime::vector<cmplx> propagation_table_;
   spacetime::vector3d<cmplx> source_table_, obs_table_;
 
+  Eigen::Array3Xcd temp_observers;
+
   TransformPair spatial_vector_transforms_;
 
   spacetime::vector<cmplx> make_propagation_table() const;
   void fill_source_table(const int);
   void propagate(const int);
+  void propagate_present_field(const int);
   void fill_results_table(const int);
 
   void fill_gmatrix_table(spacetime::vector<cmplx> &) const;