refactor: Refactor public API of collision operator

include/lbt/populations/collision/bgk.hpp

@@ -8,20 +8,53 @@
 #define LBT__COLLISION__BGK
 #pragma once
 
+#include <memory>
+
+#include "lbt/continuums/continuum_base.hpp"
+#include "lbt/populations/collision/collision_base.hpp"
+#include "lbt/populations/indexing/timestep.hpp"
+#include "lbt/units/characteristic_numbers.hpp"
+#include "lbt/converter.hpp"
+
 
 namespace lbt {
   namespace collision {
 
     /**\class  Bgk
-     * \brief  Class that holds macroscopic values with indexing based on the A-A access pattern
+     * \brief  Bhatnagar-Gross-Krook (BGK) collision operator for arbitrary lattice
      *
-     * \tparam LT   Static lattice::DdQq class containing discretisation parameters
+     * \note   "A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged
+     *          and Neutral One-Component Systems"
+     *         P.L. Bhatnagar, E.P. Gross, M. Krook
+     *         Physical Review 94 (1954)
+     *         DOI: 10.1103/PhysRev.94.511
+     *
+     * \tparam LT     Static lattice::DdQq class containing discretisation parameters
     */
+
     template <class LT>
-    class Bgk {
+    class Bgk: public CollisionBase<LT> {
       public:
+        Bgk(std::shared_ptr<lbt::Population<LT>> population, std::shared_ptr<lbt::ContinuumBase<typename LT::type>> continuum,
+          std::shared_ptr<lbt::Converter> converter, lbt::unit::ReynoldsNumber const& reynolds_number,
+          double const lbm_velocity, double const lbm_length);
+        Bgk() = default;
+        Bgk(Bgk const&) = default;
+        Bgk& operator= (Bgk const&) = default;
+        Bgk(Bgk&&) = default;
+        Bgk& operator= (Bgk&&) = default;
     };
 
+    template <class LT>
+    Bgk<LT>::Bgk(std::shared_ptr<lbt::Population<LT>> population, std::shared_ptr<lbt::ContinuumBase<typename LT::type>> continuum,
+      std::shared_ptr<lbt::Converter> converter, lbt::unit::ReynoldsNumber const& reynolds_number,
+      double const lbm_velocity, double const lbm_length)
+    : CollisionBase<LT>{population, continuum, converter} {
+      return;
+    }
+
+    // TODO(tobit): Implement collide-stream
+
   }
 
 }

include/lbt/populations/collision/collision_base.hpp

@@ -0,0 +1,77 @@
+/**
+ * \file     collision_base.hpp
+ * \mainpage Base class for all collision operators
+ * \author   Tobit Flatscher (github.com/2b-t)
+*/
+
+#ifndef LBT__COLLISION__COLLISION_BASE
+#define LBT__COLLISION__COLLISION_BASE
+#pragma once
+
+#include <memory>
+
+#include "lbt/continuums/continuum_base.hpp"
+#include "lbt/populations/indexing/timestep.hpp"
+#include "lbt/populations/population.hpp"
+#include "lbt/units/characteristic_numbers.hpp"
+#include "lbt/converter.hpp"
+
+
+namespace lbt {
+  namespace collision {
+
+    /**\class  CollisionBase
+     * \brief  Base class for all collision operators
+     *
+     * \tparam LT   Static lattice::DdQq class containing discretisation parameters
+    */
+    template <class LT>
+    class CollisionBase {
+      public:
+        template<lbt::Timestep TS>
+        void initialize();
+
+        template<lbt::Timestep TS>
+        void collideAndStream(bool const is_save);
+
+      protected:
+        // Allow population and unit converter to be nullptrs
+        CollisionBase(std::shared_ptr<lbt::Population<LT>> population,
+          std::shared_ptr<lbt::ContinuumBase<typename LT::type>> continuum,
+          std::shared_ptr<lbt::Converter> converter);
+        CollisionBase() = delete;
+        CollisionBase(CollisionBase const&) = default;
+        CollisionBase& operator= (CollisionBase const&) = default;
+        CollisionBase(CollisionBase&&) = default;
+        CollisionBase& operator= (CollisionBase&&) = default;
+
+        std::shared_ptr<lbt::Population<LT>> population_;
+        std::shared_ptr<lbt::ContinuumBase<typename LT::type>> continuum_;
+        std::shared_ptr<lbt::Converter> converter_;
+    };
+
+    template <class LT>
+    CollisionBase<LT>::CollisionBase(std::shared_ptr<lbt::Population<LT>> population,
+      std::shared_ptr<lbt::ContinuumBase<typename LT::type>> continuum,
+      std::shared_ptr<lbt::Converter> converter)
+    : population_{population}, continuum_{continuum}, converter_{converter} {
+      return;
+    }
+
+    template <class LT> template<lbt::Timestep TS>
+    void CollisionBase<LT>::initialize() {
+      // TODO(tobit): Initialize population
+      return;
+    }
+
+    template <class LT> template<lbt::Timestep TS>
+    void CollisionBase<LT>::collideAndStream(bool const is_save) {
+      // Implement collide stream through CRTP
+      return;
+    }
+
+  }
+
+}
+
+#endif // LBT__COLLISION__COLLISION_BASE

src/simple_main.cpp

@@ -25,11 +25,11 @@ int main(int argc, char* argv[]) {
   using namespace lbt::literals;
 
   // Flow and material properties
-  constexpr auto velocity = 1.7_mps;
-  constexpr auto length = 1.2_m;
+  constexpr auto velocity {1.7_mps};
+  constexpr auto length {1.2_m};
   using Material = lbt::material::Air;
-  constexpr auto temperature = 20.0_deg;
-  constexpr auto pressure = 1.0_atm;
+  constexpr auto temperature {20.0_deg};
+  constexpr auto pressure {1.0_atm};
   constexpr auto density {Material::equationOfState(temperature, pressure)};
   constexpr auto kinematic_viscosity {Material::kinematicViscosity(temperature, pressure)};
   constexpr lbt::unit::ReynoldsNumber reynolds_number {velocity, length, kinematic_viscosity};
@@ -40,54 +40,60 @@ int main(int argc, char* argv[]) {
   using T = double;
   using Lattice = lbt::lattice::D2Q9<T>;
   using CollisionOperator = lbt::collision::Bgk<Lattice>;
-  constexpr auto lbm_speed_of_sound {lbt::lattice::D2Q9<T>::CS};
-  constexpr auto lbm_velocity = 0.5*lbm_speed_of_sound;
-  constexpr long double lbm_density {1.0};
+  constexpr auto lbm_speed_of_sound {Lattice::CS};
+  constexpr auto lbm_velocity {0.5*lbm_speed_of_sound};
+  constexpr double lbm_density {1.0};
   constexpr std::int32_t NX {100};
   constexpr std::int32_t NY {100};
   constexpr std::int32_t NZ {100};
+  constexpr double lbm_length {static_cast<long double>(NX)};
   auto const output_path {std::filesystem::current_path()};
 
-  constexpr lbt::Converter unit_converter {
-    length, static_cast<long double>(NX),
+  auto const unit_converter {std::make_shared<lbt::Converter>(
+    length, lbm_length,
     velocity, lbm_velocity,
     density, lbm_density
-  };
+  )};
+
+  // Import geometry either from scenario or using VTK
 
   auto& omp_manager = lbt::OpenMpManager::getInstance();
   omp_manager.setThreadsNum(lbt::OpenMpManager::getThreadsMax());
 
-  // Geometry and boundary conditions
-  // std::vector<BoundaryCondition> boundary_conditions {};
-  // Append boundary conditions, allow specifying range for boundary condition easily
-
   // Continuum, population and collision operator
   auto continuum {std::make_shared<lbt::Continuum<T>>(NX, NY, NZ, output_path)};
+  continuum->initializeUniform(pressure, velocity, 0.0_mps, 0.0_mps);
+  // continuum->initialize(lambda)
+  // Initialize with lambda, allow user to use lambda with normalized coordinates [0, 1]
+  // Prepare lambdas for uniform and Poiseuille profile etc.
+
   auto population {std::make_shared<lbt::Population<Lattice>>(NX, NY, NZ)};
-  // auto CollisionOperator collision_operator {continuum, population, reynolds_number, lbm_velocity, lbm_length};
-  // Writer class for writing results to file
+  auto collision_operator {std::make_shared<CollisionOperator>(
+    population, continuum, unit_converter,
+    reynolds_number, lbm_velocity, lbm_length
+  )};
+  collision_operator->initialize<lbt::Timestep::Even>();
 
-  // collision_operator.initializeUniform<lbt::Timestep::Even>({u, v, w}, lbm_density);
-  // collision_operator.initializeFromContinuum<lbt::Timestep::Even>(continuum);
+  // Geometry and boundary conditions
+  // std::vector<BoundaryCondition> boundary_conditions {};
+  // Append boundary conditions, allow specifying range for boundary condition easily
 
   // Convert time-steps to simulation time
-  constexpr auto NT {static_cast<std::int64_t>(unit_converter.toLbm(simulation_time))};
-  constexpr auto NT_SAVE {static_cast<std::int64_t>(unit_converter.toLbm(save_time_step))};
+  auto const NT {static_cast<std::int64_t>(unit_converter->toLbm(simulation_time))};
+  auto const NT_SAVE {static_cast<std::int64_t>(unit_converter->toLbm(save_time_step))};
   for (std::int64_t i = 0; i < NT; i += 2) {
     // Loop performs two combined iterations at once
 
-    // TODO(tobit): Currently saving every 100th step
-    bool const is_save {(i % NT_SAVE) == 0};
-
-    // collision_operator.collideAndStream<lbt::Timestep::Even>(is_save);
+    bool is_save {(i % NT_SAVE) == 0};
+    collision_operator->collideAndStream<lbt::Timestep::Even>(is_save);
     // Enforce boundary condition on continuum
 
-    // collision_operator.collideAndStream<lbt::Timestep::Odd>(is_save);
+    is_save = (((i + 1) % NT_SAVE) == 0);
+    collision_operator->collideAndStream<lbt::Timestep::Odd>(is_save);
     // Enforce boundary condition on continuum
 
     if (is_save) {
-      auto const simulation_time {unit_converter.toPhysical<lbt::unit::Time>(i)};
-      // Convert units from LBM units
+      auto const simulation_time {unit_converter->toPhysical<lbt::unit::Time>(i)};
       continuum->save(static_cast<double>(simulation_time.get()));
     }
   }