Allow customization of Grid2D's physical dimensionalization

2 files changed, 99 insertions, 27 deletions

diff --git a/src/refinement/grid2D.h b/src/refinement/grid2D.h
index 95324f7..c90078b 100644
--- a/src/refinement/grid2D.h
+++ b/src/refinement/grid2D.h
@@ -66,10 +66,29 @@ using RelaxationTime = NamedType<T,struct NamedRelaxationTime>;
 template <typename T>
 using LatticeVelocity = NamedType<T,struct NamedLatticeVelocity>;
 
+template <typename T>
+struct Characteristics {
+  Characteristics(T l, T u, T nu, T rho):
+    length(l),
+    velocity(u),
+    viscosity(nu),
+    density(rho) { }
+
+  Characteristics(int Re):
+    Characteristics(1.0, 1.0, 1.0/Re, 1.0) { }
+
+  const T length;
+  const T velocity;
+  const T viscosity;
+  const T density;
+};
+
+
 template <typename T, template<typename> class DESCRIPTOR>
 class Grid2D {
 protected:
   FunctorPtr<IndicatorF2D<T>> _domainF;
+  const Characteristics<T> _characteristics;
 
   std::unique_ptr<UnitConverter<T,DESCRIPTOR>>  _converter;
   std::unique_ptr<CuboidGeometry2D<T>>          _cuboids;
@@ -79,6 +98,7 @@ protected:
 
   std::vector<std::unique_ptr<Dynamics<T,DESCRIPTOR>>> _dynamics;
   std::vector<std::unique_ptr<sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>>> _onLatticeBoundaryConditions;
+  std::vector<std::unique_ptr<sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>>> _offLatticeBoundaryConditions;
 
   std::vector<std::unique_ptr<RefiningGrid2D<T,DESCRIPTOR>>> _fineGrids;
 
@@ -86,8 +106,19 @@ protected:
   std::vector<std::unique_ptr<CoarseCoupler2D<T,DESCRIPTOR>>> _coarseCouplers;
 
 public:
-  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, int resolution, RelaxationTime<T> tau, int re);
-  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, int resolution, LatticeVelocity<T> uMax, int re);
+  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF,
+         RelaxationTime<T> tau,
+         int resolution,
+         Characteristics<T> characteristics);
+  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF,
+         LatticeVelocity<T> latticeVelocity,
+         int resolution,
+         Characteristics<T> characteristics);
+
+  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, RelaxationTime<T> tau,   int resolution, int re);
+  Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF, LatticeVelocity<T> uMax, int resolution, int re);
+
+  Characteristics<T> getCharacteristics() const;
 
   UnitConverter<T,DESCRIPTOR>& getConverter();
   CuboidGeometry2D<T>& getCuboidGeometry();
@@ -97,6 +128,7 @@ public:
 
   Dynamics<T,DESCRIPTOR>& addDynamics(std::unique_ptr<Dynamics<T,DESCRIPTOR>>&& dynamics);
   sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& getOnLatticeBoundaryCondition();
+  sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& getOffLatticeBoundaryCondition();
 
   void collideAndStream();
 
diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh
index 504808a..6e9405e 100644
--- a/src/refinement/grid2D.hh
+++ b/src/refinement/grid2D.hh
@@ -34,19 +34,19 @@ namespace olb {
 
 
 template <typename T, template<typename> class DESCRIPTOR>
-Grid2D<T,DESCRIPTOR>::Grid2D(
-  FunctorPtr<IndicatorF2D<T>>&& domainF,
-  int resolution,
-  RelaxationTime<T> tau,
-  int re):
+Grid2D<T,DESCRIPTOR>::Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF,
+                             RelaxationTime<T> tau,
+                             int resolution,
+                             Characteristics<T> characteristics):
   _domainF(std::move(domainF)),
+  _characteristics(characteristics),
   _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
                resolution, // resolution: number of voxels per charPhysL
                tau,        // latticeRelaxationTime: relaxation time, has to be greater than 0.5!
-               T{1},       // charPhysLength: reference length of simulation geometry
-               T{1},       // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
-               T{1./re},   // physViscosity: physical kinematic viscosity in __m^2 / s__
-               T{1})),     // physDensity: physical density in __kg / m^3__
+               characteristics.length,     // charPhysLength: reference length of simulation geometry
+               characteristics.velocity,   // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
+               characteristics.viscosity,  // physViscosity: physical kinematic viscosity in __m^2 / s__
+               characteristics.density)),  // physDensity: physical density in __kg / m^3__
   _cuboids(new CuboidGeometry2D<T>(
              *_domainF,
              _converter->getConversionFactorLength(),
@@ -69,19 +69,19 @@ Grid2D<T,DESCRIPTOR>::Grid2D(
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-Grid2D<T,DESCRIPTOR>::Grid2D(
-  FunctorPtr<IndicatorF2D<T>>&& domainF,
-  int resolution,
-  LatticeVelocity<T> velocity,
-  int re):
+Grid2D<T,DESCRIPTOR>::Grid2D(FunctorPtr<IndicatorF2D<T>>&& domainF,
+                             LatticeVelocity<T> latticeVelocity,
+                             int resolution,
+                             Characteristics<T> characteristics):
   _domainF(std::move(domainF)),
+  _characteristics(characteristics),
   _converter(new UnitConverterFromResolutionAndLatticeVelocity<T,DESCRIPTOR>(
-               resolution, // resolution: number of voxels per charPhysL
-               velocity,   // maxLatticeVelocity
-               T{1},       // charPhysLength: reference length of simulation geometry
-               T{1},       // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
-               T{1./re},   // physViscosity: physical kinematic viscosity in __m^2 / s__
-               T{1})),     // physDensity: physical density in __kg / m^3__
+               resolution,      // resolution: number of voxels per charPhysL
+               latticeVelocity, // charLatticeVelocity
+               characteristics.length,     // charPhysLength: reference length of simulation geometry
+               characteristics.velocity,   // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
+               characteristics.viscosity,  // physViscosity: physical kinematic viscosity in __m^2 / s__
+               characteristics.density)),  // physDensity: physical density in __kg / m^3__
   _cuboids(new CuboidGeometry2D<T>(
              *_domainF,
              _converter->getConversionFactorLength(),
@@ -104,6 +104,36 @@ Grid2D<T,DESCRIPTOR>::Grid2D(
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
+Grid2D<T,DESCRIPTOR>::Grid2D(
+  FunctorPtr<IndicatorF2D<T>>&& domainF,
+  RelaxationTime<T> tau,
+  int resolution,
+  int re):
+  Grid2D(std::forward<decltype(domainF)>(domainF),
+         tau,
+         resolution,
+         Characteristics<T>(re))
+{ }
+
+template <typename T, template<typename> class DESCRIPTOR>
+Grid2D<T,DESCRIPTOR>::Grid2D(
+  FunctorPtr<IndicatorF2D<T>>&& domainF,
+  LatticeVelocity<T> latticeVelocity,
+  int resolution,
+  int re):
+  Grid2D(std::forward<decltype(domainF)>(domainF),
+         latticeVelocity,
+         resolution,
+         Characteristics<T>(re))
+{ }
+
+template <typename T, template<typename> class DESCRIPTOR>
+Characteristics<T> Grid2D<T,DESCRIPTOR>::getCharacteristics() const
+{
+  return _characteristics;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
 UnitConverter<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getConverter()
 {
   return *_converter;
@@ -151,6 +181,14 @@ sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getOnLatticeB
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
+sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getOffLatticeBoundaryCondition()
+{
+  _offLatticeBoundaryConditions.emplace_back(
+    new sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>(getSuperLattice()));
+  return *_offLatticeBoundaryConditions.back();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
 void Grid2D<T,DESCRIPTOR>::collideAndStream()
 {
   for ( auto& fineCoupler : _fineCouplers ) {
@@ -214,7 +252,10 @@ template <typename T, template<typename> class DESCRIPTOR>
 Vector<T,2> Grid2D<T,DESCRIPTOR>::alignExtendToGrid(Vector<T,2> extend) const
 {
   const T deltaX = _converter->getPhysDeltaX();
-  return util::floor(extend / deltaX) * deltaX;
+  return {
+    static_cast<int>(std::floor(extend[0] / deltaX)) * deltaX,
+    static_cast<int>(std::floor(extend[1] / deltaX)) * deltaX
+  };
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -273,13 +314,12 @@ RefiningGrid2D<T,DESCRIPTOR>::RefiningGrid2D(
   Grid2D<T,DESCRIPTOR>& parentGrid, Vector<T,2> origin, Vector<T,2> extend):
   Grid2D<T,DESCRIPTOR>(
     std::unique_ptr<IndicatorF2D<T>>(new IndicatorCuboid2D<T>(extend, origin)),
-    2*parentGrid.getConverter().getResolution(),
     RelaxationTime<T>(2*parentGrid.getConverter().getLatticeRelaxationTime() - 0.5),
-    parentGrid.getConverter().getReynoldsNumber()),
+    2*parentGrid.getConverter().getResolution(),
+    parentGrid.getCharacteristics()),
   _origin(origin),
   _extend(extend),
-  _parentGrid(parentGrid)
-{ }
+  _parentGrid(parentGrid) { }
 
 template <typename T, template<typename> class DESCRIPTOR>
 Vector<T,2> RefiningGrid2D<T,DESCRIPTOR>::getOrigin() const