Feature update v0.2.0, Instrumentation, 1D Fourier

Feature update v.0.2.0

.travis.yml

@@ -7,7 +7,7 @@ language: cpp
 
 env:
  global:
-  - DOCKER_TAG=v9.1.1
+  - DOCKER_TAG=v9.2.0-1
 
 services:
  - docker
@@ -25,7 +25,7 @@ install:
  - sudo docker run $ci_env -it -u="root" -d -w="/home/bart/BART" --name bart jsrehak/bart:$DOCKER_TAG
 
 script:
- - docker exec bart cmake .
+ - docker exec bart cmake -DCMAKE_BUILD_TYPE=Debug .
  - docker exec bart make -j4
  - docker exec bart bash -c "./bart_test"
  - docker exec bart bash -c "mpirun -np 2 --allow-run-as-root --oversubscribe ./bart_test --mpi -l 0"

CMakeLists.txt

@@ -16,7 +16,11 @@
 CMAKE_MINIMUM_REQUIRED(VERSION 3.0)
 PROJECT(bart)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
+
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release)
+endif()
 
 ### DEPENDENCIES #####################################################
 # Check that DEAL II is installed
@@ -66,6 +70,9 @@ get_filename_component(PROTO_INC_DIR ${PROTO_HDRS} DIRECTORY)
 SET(DEAL_II_CXX_FLAGS_DEBUG "${DEAL_II_CXX_FLAGS_DEBUG} -g -O0 --coverage")
 SET(DEAL_II_LINKER_FLAGS_DEBUG "${DEAL_II_LINKER_FLAGS_DEBUG} --coverage")
 
+### OTHER DEPENDENCIES ###############################################
+find_package(fmt)
+
 ######################################################################
 
 DEAL_II_INITIALIZE_CACHED_VARIABLES()
@@ -98,10 +105,11 @@ ADD_EXECUTABLE(bart_test ${testing_sources})
 # Add testing definition and library to bart_test
 target_compile_definitions(bart_test PUBLIC -DTEST)
 
-
-target_link_libraries(bart ${Protobuf_LIBRARIES})
+target_link_libraries(bart ${Protobuf_LIBRARIES} -lfftw3 fmt::fmt)
 target_link_libraries(bart_test ${GTEST_BOTH_LIBRARIES}
-  ${GMOCK_BOTH_LIBRARIES} ${Protobuf_LIBRARIES})
+  ${GMOCK_BOTH_LIBRARIES} ${Protobuf_LIBRARIES} -lfftw3 fmt::fmt)
+
+
 
 DEAL_II_SETUP_TARGET(bart)
 DEAL_II_SETUP_TARGET(bart_test)

README.md

@@ -44,6 +44,8 @@ One of the major design goals of BART is to provide a framework for testing acce
 - Nonlinear diffusion acceleration (NDA)
 - Two-grid acceleration (TG)
 
+Calculation of discrete Fourier transforms is provided via the [FFTW subroutine library](http://www.fftw.org/).
+
 ### Benchmarks
 
 Benchmarks from Sood (1999) are provided in the `benchmarks` folder for validation of the code.

benchmarks/picca_2016/figure_2_saaf_dft.prm

@@ -0,0 +1,22 @@
+set problem dimension                        = 1
+set transport model                          = saaf
+set number of groups                         = 1
+set angular quadrature name                  = gauss_legendre
+set angular quadrature order                 = 4
+set do eigenvalue calculations               = true
+set reflective boundary names                =
+set do dft of error                          = true
+
+set x, y, z max values of boundary locations = 100.0
+set number of cells for x, y, z directions   = 20
+set uniform refinements                      = 2
+set number of materials                      = 2
+
+set finite element polynomial degree         = 1
+
+set output file name base                    = figure_2_saaf
+
+subsection material ID map
+set material id file name map                = 1: reactor, 2: reflector
+set material id file name                    = figure_2.material_map
+end

src/calculator/fourier/fourier_transform_fftw.cc

@@ -0,0 +1,48 @@
+#include "calculator/fourier/fourier_transform_fftw.h"
+
+namespace bart {
+
+namespace calculator {
+
+namespace fourier {
+
+FourierTransformFFTW::FourierTransformFFTW(const int n_samples)
+    : n_samples_(n_samples) {
+  input_.resize(n_samples);
+  output_.resize(n_samples);
+  input_ptr_ = reinterpret_cast<fftw::fftw_complex*>(input_.data());
+  output_ptr_ = reinterpret_cast<fftw::fftw_complex*>(output_.data());
+  plan_ = fftw::fftw_plan_dft_1d(n_samples_, input_ptr_, output_ptr_,
+                                 FFTW_FORWARD, FFTW_ESTIMATE_PATIENT);
+}
+
+FourierTransformFFTW::~FourierTransformFFTW() {
+  fftw::fftw_destroy_plan(plan_);
+}
+
+std::vector<std::complex<double>> FourierTransformFFTW::CalculateDFT(
+    const std::vector<std::complex<double>>& input,
+    Normalized normalized) {
+  input_ = input;
+  fftw::fftw_execute(plan_);
+  if (normalized) {
+    for (auto& value : output_)
+      value /= n_samples_;
+  }
+  return output_;
+}
+std::vector<std::complex<double>> FourierTransformFFTW::CalculateDFT(
+    const dealii::Vector<double>& input,
+    Normalized normalized) {
+  std::vector<std::complex<double>> input_copy(n_samples_);
+  for (int i = 0; i < n_samples_; ++i) {
+    input_copy.at(i).real(input[i]);
+  }
+  return CalculateDFT(input_copy, normalized);
+}
+
+} // namespace fourier
+
+} // namespace calculator
+
+} // namespace bart

src/calculator/fourier/fourier_transform_fftw.h

@@ -0,0 +1,48 @@
+#ifndef BART_SRC_CALCULATOR_FOURIER_FOURIER_TRANSFORM_FFTW_H_
+#define BART_SRC_CALCULATOR_FOURIER_FOURIER_TRANSFORM_FFTW_H_
+
+#include <complex>
+#include <vector>
+
+#include "calculator/fourier/fourier_transform_i.h"
+#include "utility/named_type.h"
+
+namespace bart {
+
+namespace calculator {
+
+namespace fourier {
+
+namespace fftw {
+#include <fftw3.h>
+} // namespace fftw
+
+class FourierTransformFFTW : public FourierTransformI {
+ public:
+
+  explicit FourierTransformFFTW(const int n_samples);
+  ~FourierTransformFFTW();
+
+  std::vector<std::complex<double>> CalculateDFT(
+      const std::vector<std::complex<double>>& input,
+      Normalized normalized = Normalized(false)) override;
+  std::vector<std::complex<double>> CalculateDFT(
+      const dealii::Vector<double> &input,
+      Normalized normalized = Normalized(false)) override;
+
+  int n_samples() const { return n_samples_; }
+ private:
+  const int n_samples_;
+  std::vector<std::complex<double>> input_, output_;
+  fftw::fftw_complex* input_ptr_;
+  fftw::fftw_complex* output_ptr_;
+  fftw::fftw_plan_s* plan_;
+};
+
+} // namespace fourier
+
+} // namespace calculator
+
+} // namespace bart
+
+#endif //BART_SRC_CALCULATOR_FOURIER_FOURIER_TRANSFORM_FFTW_H_

src/calculator/fourier/fourier_transform_i.h

@@ -0,0 +1,36 @@
+#ifndef BART_SRC_CALCULATOR_FOURIER_FOURIER_TRANSFORM_I_H_
+#define BART_SRC_CALCULATOR_FOURIER_FOURIER_TRANSFORM_I_H_
+
+#include <complex>
+#include <vector>
+
+#include <deal.II/lac/vector.h>
+
+#include "utility/named_type.h"
+
+namespace bart {
+
+namespace calculator {
+
+namespace fourier {
+
+using Normalized = utility::NamedType<bool, struct NormalizedStruct>;
+
+class FourierTransformI {
+ public:
+  virtual ~FourierTransformI() = default;
+  virtual std::vector<std::complex<double>> CalculateDFT(
+      const std::vector<std::complex<double>>& input,
+      Normalized) = 0;
+  virtual std::vector<std::complex<double>> CalculateDFT(
+      const dealii::Vector<double>& input,
+      Normalized) = 0;
+};
+
+} // namespace fourier
+
+} // namespace calculator
+
+} // namespace bart
+
+#endif //BART_SRC_CALCULATOR_FOURIER_FOURIER_TRANSFORM_I_H_

src/calculator/fourier/tests/fourier_transform_fftw_test.cc

@@ -0,0 +1,100 @@
+#include "calculator/fourier/fourier_transform_fftw.h"
+#include "test_helpers/gmock_wrapper.h"
+
+#include <cmath>
+#include <complex.h>
+
+namespace  {
+
+using namespace bart;
+using ::testing::ContainerEq;
+namespace fftw = bart::calculator::fourier::fftw;
+
+class CalculatorFourierTransformFFTWTest : public ::testing::Test {
+ public:
+  using FourierTransformType = calculator::fourier::FourierTransformFFTW;
+  void SetUp() override;
+  static constexpr int n_points{1000};
+  std::vector<std::complex<double>> function_, expected_fourier_transform_,
+      normalized_expected_fourier_transform_;
+  std::unique_ptr<FourierTransformType> test_transformer_ptr_;
+};
+
+void CalculatorFourierTransformFFTWTest::SetUp() {
+  test_transformer_ptr_ = std::make_unique<FourierTransformType>(n_points);
+  function_.resize(n_points);
+  expected_fourier_transform_.resize(n_points);
+  normalized_expected_fourier_transform_.resize(n_points);
+
+  auto forward_input = reinterpret_cast<fftw::fftw_complex*>(function_.data());
+  auto forward_output = reinterpret_cast<fftw::fftw_complex*>(expected_fourier_transform_.data());
+  auto forwards_plan = fftw::fftw_plan_dft_1d(n_points, forward_input,
+                                              forward_output, FFTW_FORWARD,
+                                              FFTW_ESTIMATE);
+  for (int i = 0; i < n_points; ++i) {
+    const double pi = M_PI;
+    const double x = 2*pi*static_cast<double>(i)/n_points;
+    function_.at(i) = std::sin(x) * std::cos(2*pi*x);
+  }
+  fftw::fftw_execute(forwards_plan);
+  fftw::fftw_destroy_plan(forwards_plan);
+
+  normalized_expected_fourier_transform_ = expected_fourier_transform_;
+  for (auto& value : normalized_expected_fourier_transform_)
+    value /= n_points;
+}
+
+TEST_F(CalculatorFourierTransformFFTWTest, Constructor) {
+  EXPECT_NO_THROW({
+    FourierTransformType test_transformer(n_points);
+                  });
+}
+
+TEST_F(CalculatorFourierTransformFFTWTest, Getters) {
+  EXPECT_EQ(test_transformer_ptr_->n_samples(), n_points);
+}
+
+TEST_F(CalculatorFourierTransformFFTWTest, CosineStdVector) {
+  auto calculated_fourier_transform =
+      test_transformer_ptr_->CalculateDFT(function_);
+  EXPECT_THAT(calculated_fourier_transform,
+              ContainerEq(expected_fourier_transform_));
+}
+
+TEST_F(CalculatorFourierTransformFFTWTest, CosineStdVectorNormalized) {
+  using bart::calculator::fourier::Normalized;
+
+  auto calculated_fourier_transform =
+      test_transformer_ptr_->CalculateDFT(function_, Normalized(true));
+
+  EXPECT_THAT(calculated_fourier_transform,
+              ContainerEq(normalized_expected_fourier_transform_));
+}
+
+TEST_F(CalculatorFourierTransformFFTWTest, CosineDealiiVector) {
+  dealii::Vector<double> function_vector(n_points);
+  for (int i = 0; i < n_points; ++i)
+    function_vector[i] = function_.at(i).real();
+
+  auto calculated_fourier_transform =
+      test_transformer_ptr_->CalculateDFT(function_vector);
+  EXPECT_THAT(calculated_fourier_transform,
+              ContainerEq(expected_fourier_transform_));
+}
+
+TEST_F(CalculatorFourierTransformFFTWTest, CosineDealiiVectorNormalized) {
+  using bart::calculator::fourier::Normalized;
+  dealii::Vector<double> function_vector(n_points);
+  for (int i = 0; i < n_points; ++i)
+    function_vector[i] = function_.at(i).real();
+
+  auto calculated_fourier_transform =
+      test_transformer_ptr_->CalculateDFT(function_vector, Normalized(true));
+
+  EXPECT_THAT(calculated_fourier_transform,
+              ContainerEq(normalized_expected_fourier_transform_));
+}
+
+
+
+} // namespace

src/calculator/fourier/tests/fourier_transform_mock.h

@@ -0,0 +1,34 @@
+#ifndef BART_SRC_CALCULATOR_FOURIER_TESTS_FOURIER_TRANSFORM_MOCK_H_
+#define BART_SRC_CALCULATOR_FOURIER_TESTS_FOURIER_TRANSFORM_MOCK_H_
+
+#include "calculator/fourier/fourier_transform_i.h"
+#include "test_helpers/gmock_wrapper.h"
+
+namespace bart {
+
+namespace calculator {
+
+namespace fourier {
+
+class FourierTransformMock : public FourierTransformI {
+ public:
+  MOCK_METHOD(std::vector<std::complex<double>>,
+              CalculateDFT,
+              (const std::vector<std::complex<double>>&, Normalized),
+              (override));
+  MOCK_METHOD(std::vector<std::complex<double>>,
+              CalculateDFT,
+              (const dealii::Vector<double>&, Normalized),
+              (override));
+};
+
+
+} // namespace fourier
+
+} // namespace calculator
+
+} // namespace bart
+
+
+
+#endif //BART_SRC_CALCULATOR_FOURIER_TESTS_FOURIER_TRANSFORM_MOCK_H_

src/convergence/final.h

@@ -4,7 +4,7 @@
 #include <optional>
 
 #include "convergence/final_i.h"
-#include "convergence/status.h"
+#include "convergence/status.hpp"
 
 namespace bart {
 

src/convergence/final_checker_or_n.cc

@@ -15,8 +15,6 @@ Status FinalCheckerOrN<CompareType, CheckerType>::CheckFinalConvergence(
     CompareType& previous_iteration) {
 
   StatusDeltaAndIterate(current_iteration, previous_iteration);
-  if (reporter_ptr_ != nullptr)
-    reporter_ptr_->Report(convergence_status_);
   return convergence_status_;
 }
 
@@ -28,8 +26,6 @@ Status FinalCheckerOrN<system::moments::MomentsMap ,
 
   StatusDeltaAndIterate(current_iteration, previous_iteration);
   convergence_status_.failed_index = checker_ptr_->failed_index();
-  if (reporter_ptr_ != nullptr)
-    reporter_ptr_->Report(convergence_status_);
   return convergence_status_;
 }