Initialize at openlb-1-3

1 files changed, 542 insertions, 0 deletions

diff --git a/src/functors/lattice/reductionF3D.hh b/src/functors/lattice/reductionF3D.hh
new file mode 100644
index 0000000..8837524
--- /dev/null
+++ b/src/functors/lattice/reductionF3D.hh
@@ -0,0 +1,542 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2012-2017 Lukas Baron, Tim Dornieden, Mathias J. Krause,
+ *  Albert Mink, Fabian Klemens, Benjamin Förster, Marie-Luise Maier,
+ *  Adrian Kummerlaender
+ *  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 REDUCTION_F_3D_HH
+#define REDUCTION_F_3D_HH
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+#include <algorithm>
+#include "functors/lattice/reductionF3D.h"
+#include "dynamics/lbHelpers.h"  // for computation of lattice rho and velocity
+
+namespace olb {
+
+template<typename T, typename DESCRIPTOR>
+SuperLatticeFfromAnalyticalF3D<T, DESCRIPTOR>::SuperLatticeFfromAnalyticalF3D(
+  FunctorPtr<AnalyticalF3D<T,T>>&& f,
+  SuperLattice3D<T, DESCRIPTOR>&   sLattice)
+  : SuperLatticeF3D<T, DESCRIPTOR>(sLattice, f->getTargetDim()),
+    _f(std::move(f))
+{
+  this->getName() = "fromAnalyticalF(" + _f->getName() + ")";
+
+  LoadBalancer<T>&     load   = sLattice.getLoadBalancer();
+  CuboidGeometry3D<T>& cuboid = sLattice.getCuboidGeometry();
+
+  for (int iC = 0; iC < load.size(); ++iC) {
+    this->_blockF.emplace_back(
+      new BlockLatticeFfromAnalyticalF3D<T,DESCRIPTOR>(
+        *_f,
+        sLattice.getExtendedBlockLattice(iC),
+        cuboid.get(load.glob(iC)))
+    );
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+bool SuperLatticeFfromAnalyticalF3D<T, DESCRIPTOR>::operator()(
+  T output[], const int input[])
+{
+  T physR[3] = {};
+  this->_sLattice.getCuboidGeometry().getPhysR(physR,input);
+  return _f(output,physR);
+}
+
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeFfromAnalyticalF3D<T, DESCRIPTOR>::BlockLatticeFfromAnalyticalF3D(
+  AnalyticalF3D<T, T>&                    f,
+  BlockLatticeStructure3D<T, DESCRIPTOR>& lattice,
+  Cuboid3D<T>&                            cuboid)
+  : BlockLatticeF3D<T, DESCRIPTOR>(lattice, f.getTargetDim()),
+    _f(f),
+    _cuboid(cuboid)
+{
+  this->getName() = "blockFfromAnalyticalF(" + _f.getName() + ")";
+}
+
+template<typename T, typename DESCRIPTOR>
+bool BlockLatticeFfromAnalyticalF3D<T, DESCRIPTOR>::operator()(
+  T output[], const int input[])
+{
+  T physR[3] = {};
+  _cuboid.getPhysR(physR,input);
+  return _f(output,physR);
+}
+
+
+//////////// not yet working // symbolically ///////////////////
+////////////////////////////////////////////////
+template<typename T, typename DESCRIPTOR>
+SmoothBlockIndicator3D<T, DESCRIPTOR>::SmoothBlockIndicator3D(
+  IndicatorF3D<T>& f, T h, T eps, int wa)
+  : BlockDataF3D<T, T>((int)((f.getMax()[0] - f.getMin()[0]) / h + (wa+1)*eps)+2,
+                       (int)((f.getMax()[1] - f.getMin()[1]) / h + (wa+1)*eps)+2,
+                       (int)((f.getMax()[2] - f.getMin()[2]) / h + (wa+1)*eps)+2),
+    _f(f),
+    _h(h),
+    _eps(eps),
+    _wa(wa)
+{
+  this->getName() = "SmoothBlockIndicator3D";
+
+  // shrink epsilon boundary to one size
+  this->_eps = this->_eps/((this->_wa-1)/2.0);
+  T value;
+
+  // Important parameter of the gaussian psf (point spread function), but adaption should not be necessary
+  T sigma = 1.0;
+
+  T weights[this->_wa][this->_wa][this->_wa];
+  T sum = 0;
+  int iStart = floor(this->_wa/2.0);
+  int iEnd = ceil(this->_wa/2.0);
+
+  // calculate weights: they are constants, but calculation here is less error-prone than hardcoding these parameters
+  for (int x = -iStart; x < iEnd; x++) {
+    for (int y = -iStart; y < iEnd; y++) {
+      for (int z = -iStart; z < iEnd; z++) {
+        weights[x+iStart][y+iStart][z+iStart] = exp(-(x*x+y*y+z*z)/(2*sigma*sigma)) / (pow(sigma,3)*sqrt(pow(2,3)*pow(M_PI,3)));
+        // important because sum of all weigths only equals 1 for this->_wa -> infinity
+        sum += weights[x+iStart][y+iStart][z+iStart];
+      }
+    }
+  }
+
+  for (int iX=0; iX<this->getBlockData().getNx(); iX++) {
+    for (int iY=0; iY<this->getBlockData().getNy(); iY++) {
+      for (int iZ=0; iZ<this->getBlockData().getNz(); iZ++) {
+        bool output[1];
+        value = 0;
+
+        // input: regarded point (centre)
+        T input[] = {
+          _f.getMin()[0] + (iX-iStart*_eps)*_h,
+          _f.getMin()[1] + (iY-iStart*_eps)*_h,
+          _f.getMin()[2] + (iZ-iStart*_eps)*_h
+        };
+
+        /*
+         * three loops to look at every point (which weight is not 0) around the regarded point
+         * sum all weighted porosities
+         */
+        for (int x = -iStart; x < iEnd; x++) {
+          for (int y = -iStart; y < iEnd; y++) {
+            for (int z = -iStart; z < iEnd; z++) {
+              // move from regarded point to point in neighborhood
+              input[0] += x*_eps*_h;
+              input[1] += y*_eps*_h;
+              input[2] += z*_eps*_h;
+
+              // get porosity
+              _f(output,input);
+
+              // sum porosity
+              value += output[0] * weights[x+iStart][y+iStart][z+iStart];
+
+              // move back to centre
+              input[0] -= x*_eps*_h;
+              input[1] -= y*_eps*_h;
+              input[2] -= z*_eps*_h;
+            }
+          }
+        }
+        /*
+         * Round to 3 decimals
+         * See above sum != 1.0, that's the reason for devision, otherwise porosity will never reach 0
+         */
+        this->getBlockData().get(iX,iY,iZ,0) = nearbyint(1000*value/sum)/1000.0;
+      }
+    }
+  }
+}
+/*
+template<typename T, typename DESCRIPTOR>
+bool SmoothBlockIndicator3D<T, DESCRIPTOR>::operator()(
+  T output[], const int input[])
+{
+  T physR[3] = {};
+  _superGeometry.getPhysR(physR,input[0],input[1],input[2] );
+  _f(output,physR);
+  return true;
+}*/
+
+
+template<typename T, typename DESCRIPTOR>
+SuperLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>::
+SuperLatticeInterpPhysVelocity3Degree3D(
+  SuperLattice3D<T, DESCRIPTOR>& sLattice, UnitConverter<T,DESCRIPTOR>& conv, int range)
+  : SuperLatticeF3D<T, DESCRIPTOR>(sLattice, 3)
+{
+  this->getName() = "Interp3DegreeVelocity";
+  int maxC = this->_sLattice.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    BlockLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>* foo =
+      new BlockLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>(
+      sLattice.getExtendedBlockLattice(iC),
+      conv,
+      &sLattice.getCuboidGeometry().get(this->_sLattice.getLoadBalancer().
+                                        glob(iC)),
+      sLattice.getOverlap(),
+      range);
+    _bLattices.push_back(foo);
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void SuperLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>::operator()(
+  T output[], const T input[], const int iC)
+{
+  _bLattices[this->_sLattice.getLoadBalancer().loc(iC)]->operator()(output,
+      input);
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>::
+BlockLatticeInterpPhysVelocity3Degree3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice, UnitConverter<T,DESCRIPTOR>& conv,
+  Cuboid3D<T>* c, int overlap, int range)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3),
+    _conv(conv),
+    _cuboid(c),
+    _overlap(overlap),
+    _range(range)
+{
+  this->getName() = "BlockLatticeInterpVelocity3Degree3D";
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>::
+BlockLatticeInterpPhysVelocity3Degree3D(
+  const BlockLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>& rhs) :
+  BlockLatticeF3D<T, DESCRIPTOR>(rhs._blockLattice, 3),
+  _conv(rhs._conv),
+  _cuboid(rhs._cuboid),
+  _overlap(rhs._overlap),
+  _range(rhs._range)
+{
+}
+
+template<typename T, typename DESCRIPTOR>
+void BlockLatticeInterpPhysVelocity3Degree3D<T, DESCRIPTOR>::operator()(
+  T output[3], const T input[3])
+{
+  T u[3], rho, volume;
+  int latIntPos[3] = {0};
+  T latPhysPos[3] = {T()};
+  _cuboid->getFloorLatticeR(latIntPos, &input[0]);
+  _cuboid->getPhysR(latPhysPos, latIntPos);
+
+  latIntPos[0] += _overlap;
+  latIntPos[1] += _overlap;
+  latIntPos[2] += _overlap;
+
+  volume=T(1);
+  for (int i = -_range; i <= _range+1; ++i) {
+    for (int j = -_range; j <= _range+1; ++j) {
+      for (int k = -_range; k <= _range+1; ++k) {
+
+        this->_blockLattice.get(latIntPos[0]+i, latIntPos[1]+j,
+                                latIntPos[2]+k).computeRhoU(rho, u);
+        for (int l = -_range; l <= _range+1; ++l) {
+          if (l != i) {
+            volume *= (input[0] - (latPhysPos[0]+ l *_cuboid->getDeltaR()))
+                      / (latPhysPos[0] + i *_cuboid->getDeltaR()
+                         - (latPhysPos[0] + l *_cuboid->getDeltaR()));
+          }
+        }
+        for (int m = -_range; m <= _range+1; ++m) {
+          if (m != j) {
+            volume *= (input[1]
+                       - (latPhysPos[1] + m *_cuboid->getDeltaR()))
+                      / (latPhysPos[1] + j * _cuboid->getDeltaR()
+                         - (latPhysPos[1] + m * _cuboid->getDeltaR()));
+          }
+        }
+        for (int n = -_range; n <= _range+1; ++n) {
+          if (n != k) {
+            volume *= (input[2]
+                       - (latPhysPos[2] + n * _cuboid->getDeltaR()))
+                      / (latPhysPos[2] + k * _cuboid->getDeltaR()
+                         - (latPhysPos[2] + n * _cuboid->getDeltaR()));
+          }
+        }
+        output[0] += u[0] * volume;
+        output[1] += u[1] * volume;
+        output[2] += u[2] * volume;
+        volume=T(1);
+      }
+    }
+  }
+
+  output[0] = _conv.getPhysVelocity(output[0]);
+  output[1] = _conv.getPhysVelocity(output[1]);
+  output[2] = _conv.getPhysVelocity(output[2]);
+}
+
+
+template<typename T, typename DESCRIPTOR>
+SuperLatticeInterpDensity3Degree3D<T, DESCRIPTOR>::SuperLatticeInterpDensity3Degree3D(
+  SuperLattice3D<T, DESCRIPTOR>& sLattice, SuperGeometry3D<T>& sGeometry,
+  UnitConverter<T,DESCRIPTOR>& conv, int range) :
+  SuperLatticeF3D<T, DESCRIPTOR>(sLattice, 3)
+{
+  this->getName() = "Interp3DegreeDensity";
+  int maxC = this->_sLattice.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int lociC = 0; lociC < maxC; lociC++) {
+    int globiC = this->_sLattice.getLoadBalancer().glob(lociC);
+
+    BlockLatticeInterpDensity3Degree3D<T, DESCRIPTOR>* foo =
+      new BlockLatticeInterpDensity3Degree3D<T, DESCRIPTOR>(
+      sLattice.getExtendedBlockLattice(lociC),
+      sGeometry.getExtendedBlockGeometry(lociC),
+      conv,
+      &sLattice.getCuboidGeometry().get(globiC),
+      sLattice.getOverlap(), range);
+    _bLattices.push_back(foo);
+
+    if (sLattice.getOverlap() <= range + 1)
+      std::cout << "lattice overlap has to be larger than (range + 1)"
+                << std::endl;
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+SuperLatticeInterpDensity3Degree3D<T, DESCRIPTOR>::~SuperLatticeInterpDensity3Degree3D()
+{
+  // first deconstruct vector elements
+  for ( auto it : _bLattices) {
+    delete it;
+  }
+  // then delete std::vector
+  _bLattices.clear();
+}
+
+template<typename T, typename DESCRIPTOR>
+void SuperLatticeInterpDensity3Degree3D<T, DESCRIPTOR>::operator()(T output[],
+    const T input[], const int globiC)
+{
+  if (this->_sLattice.getLoadBalancer().rank(globiC) == singleton::mpi().getRank()) {
+    _bLattices[this->_sLattice.getLoadBalancer().loc(globiC)]->operator()(output,
+        input);
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeInterpDensity3Degree3D<T, DESCRIPTOR>::BlockLatticeInterpDensity3Degree3D(
+  BlockLatticeStructure3D<T, DESCRIPTOR>& blockLattice,
+  BlockGeometryStructure3D<T>& blockGeometry, UnitConverter<T,DESCRIPTOR>& conv,
+  Cuboid3D<T>* c, int overlap, int range) :
+  BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3), _blockGeometry(blockGeometry),
+  _conv(conv), _cuboid(c), _overlap(overlap), _range(range)
+{
+  this->getName() = "BlockLatticeInterpDensity3Degree3D";
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeInterpDensity3Degree3D<T, DESCRIPTOR>::BlockLatticeInterpDensity3Degree3D(
+  const BlockLatticeInterpDensity3Degree3D<T, DESCRIPTOR>& rhs) :
+  BlockLatticeF3D<T, DESCRIPTOR>(rhs._blockLattice, 3),
+  _blockGeometry(rhs._blockGeometry),_conv(rhs._conv), _cuboid(
+    rhs._cuboid), _overlap(rhs._overlap), _range(rhs._range)
+{
+}
+
+template<typename T, typename DESCRIPTOR>
+void BlockLatticeInterpDensity3Degree3D<T, DESCRIPTOR>::operator()(
+  T output[DESCRIPTOR::q], const T input[3])
+{
+  T volume = T(1);
+  T f_iPop = 0.;
+  /** neighbor position on grid of input value in lattice units
+   *referred to local cuboid
+   */
+  int latIntPos[3] = { 0 };
+  // neighbor position on grid of input value in physical units
+  T latPhysPos[3] = { T() };
+  // input is physical position on grid
+  _cuboid->getFloorLatticeR(latIntPos, input);
+  // latPhysPos is global physical position on geometry
+  _cuboid->getPhysR(latPhysPos, latIntPos);
+
+  // point on cuboid equals cell on blocklattice (extended) shifted by _overlap
+  latIntPos[0] += _overlap;
+  latIntPos[1] += _overlap;
+  latIntPos[2] += _overlap;
+
+  for (unsigned iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
+    output[iPop] = T(0);
+    for (int i = -_range; i <= _range + 1; ++i) {
+      for (int j = -_range; j <= _range + 1; ++j) {
+        for (int k = -_range; k <= _range + 1; ++k) {
+          f_iPop = 0.;
+          // just if material of cell != 1 there may be information of fluid density
+          if (_blockGeometry.getMaterial(latIntPos[0] + i, latIntPos[1] + j,
+                                         latIntPos[2] + k) != 0) {
+            // because of communication it is possible to get density information
+            // from neighboring cuboid
+            f_iPop = this->_blockLattice.get(latIntPos[0] + i, latIntPos[1] + j,
+                                             latIntPos[2] + k)[iPop];
+          }
+          for (int l = -_range; l <= _range + 1; ++l) {
+            if (l != i) {
+              volume *= (input[0] - (latPhysPos[0] + l * _cuboid->getDeltaR()))
+                        / (latPhysPos[0] + i * _cuboid->getDeltaR()
+                           - (latPhysPos[0] + l * _cuboid->getDeltaR()));
+            }
+          }
+          for (int m = -_range; m <= _range + 1; ++m) {
+            if (m != j) {
+              volume *= (input[1] - (latPhysPos[1] + m * _cuboid->getDeltaR()))
+                        / (latPhysPos[1] + j * _cuboid->getDeltaR()
+                           - (latPhysPos[1] + m * _cuboid->getDeltaR()));
+            }
+          }
+          for (int n = -_range; n <= _range + 1; ++n) {
+            if (n != k) {
+              volume *= (input[2] - (latPhysPos[2] + n * _cuboid->getDeltaR()))
+                        / (latPhysPos[2] + k * _cuboid->getDeltaR()
+                           - (latPhysPos[2] + n * _cuboid->getDeltaR()));
+            }
+          }
+          output[iPop] += f_iPop * volume;
+          volume = T(1);
+        }
+      }
+    }
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+SuperLatticeSmoothDiracDelta3D<T, DESCRIPTOR>::SuperLatticeSmoothDiracDelta3D(
+  SuperLattice3D<T, DESCRIPTOR>& sLattice,
+  UnitConverter<T,DESCRIPTOR>& conv, SuperGeometry3D<T>& sGeometry) :
+  SuperLatticeF3D<T, DESCRIPTOR>(sLattice, 3)
+{
+  this->getName() = "SuperLatticeSmoothDiracDelta3D";
+  int maxC = this->_sLattice.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int lociC = 0; lociC < maxC; lociC++) {
+    int globiC = this->_sLattice.getLoadBalancer().glob(lociC);
+
+    BlockLatticeSmoothDiracDelta3D<T, DESCRIPTOR>* foo =
+      new BlockLatticeSmoothDiracDelta3D<T, DESCRIPTOR>(
+      sLattice.getExtendedBlockLattice(lociC),
+      conv, &sLattice.getCuboidGeometry().get(globiC)
+    );
+    _bLattices.push_back(foo);
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+SuperLatticeSmoothDiracDelta3D<T, DESCRIPTOR>::~SuperLatticeSmoothDiracDelta3D()
+{
+  for ( auto it : _bLattices) {
+    delete it;
+  }
+  _bLattices.clear();
+}
+
+template<typename T, typename DESCRIPTOR>
+void SuperLatticeSmoothDiracDelta3D<T, DESCRIPTOR>::operator()(T delta[4][4][4],
+    const T physPos[3], const int globiC)
+{
+  if (this->_sLattice.getLoadBalancer().rank(globiC) == singleton::mpi().getRank()) {
+    _bLattices[this->_sLattice.getLoadBalancer().loc(globiC)]->operator()(delta,
+        physPos);
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeSmoothDiracDelta3D<T, DESCRIPTOR>::BlockLatticeSmoothDiracDelta3D(
+  BlockLattice3D<T, DESCRIPTOR>& blockLattice, UnitConverter<T,DESCRIPTOR>& conv, Cuboid3D<T>* cuboid)
+  : BlockLatticeF3D<T, DESCRIPTOR>(blockLattice, 3), _conv(conv), _cuboid(cuboid)
+{
+  this->getName() = "BlockLatticeSmoothDiracDelta3D";
+}
+
+template<typename T, typename DESCRIPTOR>
+BlockLatticeSmoothDiracDelta3D<T, DESCRIPTOR>::BlockLatticeSmoothDiracDelta3D(
+  const BlockLatticeSmoothDiracDelta3D<T, DESCRIPTOR>& rhs)
+  :
+  BlockLatticeF3D<T, DESCRIPTOR>(rhs._blockLattice, 3), _conv(rhs._conv), _cuboid(rhs._cuboid)
+{
+}
+
+template<typename T, typename DESCRIPTOR>
+void BlockLatticeSmoothDiracDelta3D<T, DESCRIPTOR>::operator()(
+  T delta[4][4][4], const T physPos[])
+{
+  int range = 1;
+  T a, b, c = T();
+  int latticeRoundedPosP[3] = { 0 };
+  T physRoundedPosP[3] = { T() };
+  T physLatticeL = _conv.getConversionFactorLength();
+
+  T counter = 0.;
+
+  _cuboid->getLatticeR(latticeRoundedPosP, physPos);
+  _cuboid->getPhysR(physRoundedPosP, latticeRoundedPosP);
+
+  for (int i = -range; i <= range + 1; ++i) {
+    for (int j = -range; j <= range + 1; ++j) {
+      for (int k = -range; k <= range + 1; ++k) {
+        delta[i+range][j+range][k+range] = T(1);
+        // a, b, c in lattice units cause physical ones get cancelled
+        a = (physRoundedPosP[0] + i * physLatticeL - physPos[0])
+            / physLatticeL;
+        b =  (physRoundedPosP[1] + j * physLatticeL - physPos[1])
+             / physLatticeL;
+        c = (physRoundedPosP[2] + k * physLatticeL - physPos[2])
+            / physLatticeL;
+
+        // the for loops already define that a, b, c are smaller than 2
+        delta[i+range][j+range][k+range] *= 1. / 4 * (1 + cos(M_PI * a / 2.));
+        delta[i+range][j+range][k+range] *= 1. / 4 * (1 + cos(M_PI * b / 2.));
+        delta[i+range][j+range][k+range] *= 1. / 4 * (1 + cos(M_PI * c / 2.));
+
+        counter += delta[i+range][j+range][k+range];
+      }
+    }
+  }
+
+  //  if (!util::nearZero(counter - T(1))){
+  //    // sum of delta has to be one
+  //    std::cout << "[" << this->getName() << "] " <<
+  //        "Delta summed up does not equal 1 but = " <<
+  //        counter << std::endl;
+  //  }
+
+}
+
+
+}  // end namespace olb
+
+#endif