Extract refinement scaffolding into separate units

1 files changed, 208 insertions, 0 deletions

diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh
new file mode 100644
index 0000000..86c3102
--- /dev/null
+++ b/src/refinement/grid2D.hh
@@ -0,0 +1,208 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2019 Adrian Kummerländer
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program is free software; you can redistribute it and/or
+ *  modify it under the terms of the GNU General Public License
+ *  as published by the Free Software Foundation; either version 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public
+ *  License along with this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+*/
+
+#ifndef REFINEMENT_GRID_2D_HH
+#define REFINEMENT_GRID_2D_HH
+
+#include "grid2D.h"
+#include "coupler2D.hh"
+
+#include "communication/heuristicLoadBalancer.h"
+
+namespace olb {
+
+
+template <typename T, template<typename> class DESCRIPTOR>
+std::unique_ptr<Grid2D<T,DESCRIPTOR>> Grid2D<T,DESCRIPTOR>::make(
+                                     IndicatorF2D<T>& domainF,
+                                     int resolution, T tau, int re)
+{
+  return std::unique_ptr<Grid2D<T,DESCRIPTOR>>(
+           new Grid2D<T,DESCRIPTOR>(domainF, resolution, tau, re)
+         );
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+Grid2D<T,DESCRIPTOR>::Grid2D(IndicatorF2D<T>& domainF, int resolution, T tau, int re):
+  _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__
+               T{1})),
+  _cuboids(new CuboidGeometry2D<T>(
+             domainF,
+             _converter->getConversionFactorLength(),
+             1)),
+  _balancer(new HeuristicLoadBalancer<T>(
+              *_cuboids)),
+  _geometry(new SuperGeometry2D<T>(
+              *_cuboids,
+              *_balancer,
+              2)),
+  _lattice(new SuperLattice2D<T,DESCRIPTOR>(
+             *_geometry))
+{
+  _converter->print();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+UnitConverter<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getConverter()
+{
+  return *_converter;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+CuboidGeometry2D<T>& Grid2D<T,DESCRIPTOR>::getCuboidGeometry()
+{
+  return *_cuboids;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+LoadBalancer<T>& Grid2D<T,DESCRIPTOR>::getLoadBalancer()
+{
+  return *_balancer;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+SuperGeometry2D<T>& Grid2D<T,DESCRIPTOR>::getSuperGeometry()
+{
+  return *_geometry;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+SuperLattice2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getSuperLattice()
+{
+  return *_lattice;
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+T Grid2D<T,DESCRIPTOR>::getScalingFactor()
+{
+  return 4. - _converter->getLatticeRelaxationFrequency();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+T Grid2D<T,DESCRIPTOR>::getInvScalingFactor()
+{
+  return 1./getScalingFactor();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+void Grid2D<T,DESCRIPTOR>::collideAndStream()
+{
+  for ( auto& fineCoupler : _fineCouplers ) {
+    fineCoupler->store();
+  }
+
+  this->getSuperLattice().collideAndStream();
+
+  for ( auto& fineGrid : _fineGrids ) {
+    fineGrid->getSuperLattice().collideAndStream();
+  }
+
+  for ( auto& fineCoupler : _fineCouplers ) {
+    fineCoupler->interpolate();
+    fineCoupler->couple();
+  }
+
+  for ( auto& fineGrid : _fineGrids ) {
+    fineGrid->getSuperLattice().collideAndStream();
+  }
+
+  for ( auto& fineCoupler : _fineCouplers ) {
+    fineCoupler->store();
+    fineCoupler->couple();
+  }
+
+  for ( auto& coarseCoupler : _coarseCouplers ) {
+    coarseCoupler->couple();
+  }
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+FineCoupler2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::addFineCoupling(
+  Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend)
+{
+  _fineCouplers.emplace_back(
+    new FineCoupler2D<T,DESCRIPTOR>(
+      *this, fineGrid, origin, extend));
+  return *_fineCouplers.back();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+CoarseCoupler2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::addCoarseCoupling(
+  Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend)
+{
+  _coarseCouplers.emplace_back(
+    new CoarseCoupler2D<T,DESCRIPTOR>(
+      *this, fineGrid, origin, extend));
+  return *_coarseCouplers.back();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(IndicatorF2D<T>& domainF)
+{
+  _fineGrids.emplace_back(
+    new Grid2D<T,DESCRIPTOR>(
+      domainF,
+      2*getConverter().getResolution(),
+      2.0*getConverter().getLatticeRelaxationTime() - 0.5,
+      getConverter().getReynoldsNumber()
+    ));
+  return *_fineGrids.back();
+}
+
+template <typename T, template<typename> class DESCRIPTOR>
+Grid2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::refine(Vector<T,2> origin, Vector<T,2> extend)
+{
+  IndicatorCuboid2D<T> fineCuboid(extend, origin);
+  auto& fineGrid = refine(fineCuboid);
+
+  const Vector<T,2> extendX = {extend[0],0};
+  const Vector<T,2> extendY = {0,extend[1]};
+
+  addFineCoupling(fineGrid, origin,           extendY);
+  addFineCoupling(fineGrid, origin + extendX, extendY);
+  addFineCoupling(fineGrid, origin + extendY, extendX);
+  addFineCoupling(fineGrid, origin,           extendX);
+
+  const T coarseDeltaX = getConverter().getPhysDeltaX();
+  const Vector<T,2> innerOrigin = origin + coarseDeltaX;
+  const Vector<T,2> innerExtendX = extendX - Vector<T,2> {2*coarseDeltaX,0};
+  const Vector<T,2> innerExtendY = extendY - Vector<T,2> {0,2*coarseDeltaX};
+
+  addCoarseCoupling(fineGrid, innerOrigin,                innerExtendY);
+  addCoarseCoupling(fineGrid, innerOrigin + innerExtendX, innerExtendY);
+  addCoarseCoupling(fineGrid, innerOrigin + innerExtendY, innerExtendX);
+  addCoarseCoupling(fineGrid, innerOrigin,                innerExtendX);
+
+  return fineGrid;
+}
+
+
+}
+
+#endif