From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/functors/lattice/reductionF3D.hh | 542 +++++++++++++++++++++++++++++++++++
1 file changed, 542 insertions(+)
create mode 100644 src/functors/lattice/reductionF3D.hh
(limited to 'src/functors/lattice/reductionF3D.hh')
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
+ *
+ *
+ * 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
+#include "functors/lattice/reductionF3D.h"
+#include "dynamics/lbHelpers.h" // for computation of lattice rho and velocity
+
+namespace olb {
+
+template
+SuperLatticeFfromAnalyticalF3D::SuperLatticeFfromAnalyticalF3D(
+ FunctorPtr>&& f,
+ SuperLattice3D& sLattice)
+ : SuperLatticeF3D(sLattice, f->getTargetDim()),
+ _f(std::move(f))
+{
+ this->getName() = "fromAnalyticalF(" + _f->getName() + ")";
+
+ LoadBalancer& load = sLattice.getLoadBalancer();
+ CuboidGeometry3D& cuboid = sLattice.getCuboidGeometry();
+
+ for (int iC = 0; iC < load.size(); ++iC) {
+ this->_blockF.emplace_back(
+ new BlockLatticeFfromAnalyticalF3D(
+ *_f,
+ sLattice.getExtendedBlockLattice(iC),
+ cuboid.get(load.glob(iC)))
+ );
+ }
+}
+
+template
+bool SuperLatticeFfromAnalyticalF3D::operator()(
+ T output[], const int input[])
+{
+ T physR[3] = {};
+ this->_sLattice.getCuboidGeometry().getPhysR(physR,input);
+ return _f(output,physR);
+}
+
+
+template
+BlockLatticeFfromAnalyticalF3D::BlockLatticeFfromAnalyticalF3D(
+ AnalyticalF3D& f,
+ BlockLatticeStructure3D& lattice,
+ Cuboid3D& cuboid)
+ : BlockLatticeF3D(lattice, f.getTargetDim()),
+ _f(f),
+ _cuboid(cuboid)
+{
+ this->getName() = "blockFfromAnalyticalF(" + _f.getName() + ")";
+}
+
+template
+bool BlockLatticeFfromAnalyticalF3D::operator()(
+ T output[], const int input[])
+{
+ T physR[3] = {};
+ _cuboid.getPhysR(physR,input);
+ return _f(output,physR);
+}
+
+
+//////////// not yet working // symbolically ///////////////////
+////////////////////////////////////////////////
+template
+SmoothBlockIndicator3D::SmoothBlockIndicator3D(
+ IndicatorF3D& f, T h, T eps, int wa)
+ : BlockDataF3D((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; iXgetBlockData().getNx(); iX++) {
+ for (int iY=0; iYgetBlockData().getNy(); iY++) {
+ for (int iZ=0; iZgetBlockData().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
+bool SmoothBlockIndicator3D::operator()(
+ T output[], const int input[])
+{
+ T physR[3] = {};
+ _superGeometry.getPhysR(physR,input[0],input[1],input[2] );
+ _f(output,physR);
+ return true;
+}*/
+
+
+template
+SuperLatticeInterpPhysVelocity3Degree3D::
+SuperLatticeInterpPhysVelocity3Degree3D(
+ SuperLattice3D& sLattice, UnitConverter& conv, int range)
+ : SuperLatticeF3D(sLattice, 3)
+{
+ this->getName() = "Interp3DegreeVelocity";
+ int maxC = this->_sLattice.getLoadBalancer().size();
+ this->_blockF.reserve(maxC);
+ for (int iC = 0; iC < maxC; iC++) {
+ BlockLatticeInterpPhysVelocity3Degree3D* foo =
+ new BlockLatticeInterpPhysVelocity3Degree3D(
+ sLattice.getExtendedBlockLattice(iC),
+ conv,
+ &sLattice.getCuboidGeometry().get(this->_sLattice.getLoadBalancer().
+ glob(iC)),
+ sLattice.getOverlap(),
+ range);
+ _bLattices.push_back(foo);
+ }
+}
+
+template
+void SuperLatticeInterpPhysVelocity3Degree3D::operator()(
+ T output[], const T input[], const int iC)
+{
+ _bLattices[this->_sLattice.getLoadBalancer().loc(iC)]->operator()(output,
+ input);
+}
+
+template
+BlockLatticeInterpPhysVelocity3Degree3D::
+BlockLatticeInterpPhysVelocity3Degree3D(
+ BlockLatticeStructure3D& blockLattice, UnitConverter& conv,
+ Cuboid3D* c, int overlap, int range)
+ : BlockLatticeF3D(blockLattice, 3),
+ _conv(conv),
+ _cuboid(c),
+ _overlap(overlap),
+ _range(range)
+{
+ this->getName() = "BlockLatticeInterpVelocity3Degree3D";
+}
+
+template
+BlockLatticeInterpPhysVelocity3Degree3D::
+BlockLatticeInterpPhysVelocity3Degree3D(
+ const BlockLatticeInterpPhysVelocity3Degree3D& rhs) :
+ BlockLatticeF3D(rhs._blockLattice, 3),
+ _conv(rhs._conv),
+ _cuboid(rhs._cuboid),
+ _overlap(rhs._overlap),
+ _range(rhs._range)
+{
+}
+
+template
+void BlockLatticeInterpPhysVelocity3Degree3D::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
+SuperLatticeInterpDensity3Degree3D::SuperLatticeInterpDensity3Degree3D(
+ SuperLattice3D& sLattice, SuperGeometry3D& sGeometry,
+ UnitConverter& conv, int range) :
+ SuperLatticeF3D(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* foo =
+ new BlockLatticeInterpDensity3Degree3D(
+ 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
+SuperLatticeInterpDensity3Degree3D::~SuperLatticeInterpDensity3Degree3D()
+{
+ // first deconstruct vector elements
+ for ( auto it : _bLattices) {
+ delete it;
+ }
+ // then delete std::vector
+ _bLattices.clear();
+}
+
+template
+void SuperLatticeInterpDensity3Degree3D::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
+BlockLatticeInterpDensity3Degree3D::BlockLatticeInterpDensity3Degree3D(
+ BlockLatticeStructure3D& blockLattice,
+ BlockGeometryStructure3D& blockGeometry, UnitConverter& conv,
+ Cuboid3D* c, int overlap, int range) :
+ BlockLatticeF3D(blockLattice, 3), _blockGeometry(blockGeometry),
+ _conv(conv), _cuboid(c), _overlap(overlap), _range(range)
+{
+ this->getName() = "BlockLatticeInterpDensity3Degree3D";
+}
+
+template
+BlockLatticeInterpDensity3Degree3D::BlockLatticeInterpDensity3Degree3D(
+ const BlockLatticeInterpDensity3Degree3D& rhs) :
+ BlockLatticeF3D(rhs._blockLattice, 3),
+ _blockGeometry(rhs._blockGeometry),_conv(rhs._conv), _cuboid(
+ rhs._cuboid), _overlap(rhs._overlap), _range(rhs._range)
+{
+}
+
+template
+void BlockLatticeInterpDensity3Degree3D::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
+SuperLatticeSmoothDiracDelta3D::SuperLatticeSmoothDiracDelta3D(
+ SuperLattice3D& sLattice,
+ UnitConverter& conv, SuperGeometry3D& sGeometry) :
+ SuperLatticeF3D(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* foo =
+ new BlockLatticeSmoothDiracDelta3D(
+ sLattice.getExtendedBlockLattice(lociC),
+ conv, &sLattice.getCuboidGeometry().get(globiC)
+ );
+ _bLattices.push_back(foo);
+ }
+}
+
+template
+SuperLatticeSmoothDiracDelta3D::~SuperLatticeSmoothDiracDelta3D()
+{
+ for ( auto it : _bLattices) {
+ delete it;
+ }
+ _bLattices.clear();
+}
+
+template
+void SuperLatticeSmoothDiracDelta3D::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
+BlockLatticeSmoothDiracDelta3D::BlockLatticeSmoothDiracDelta3D(
+ BlockLattice3D& blockLattice, UnitConverter& conv, Cuboid3D* cuboid)
+ : BlockLatticeF3D(blockLattice, 3), _conv(conv), _cuboid(cuboid)
+{
+ this->getName() = "BlockLatticeSmoothDiracDelta3D";
+}
+
+template
+BlockLatticeSmoothDiracDelta3D::BlockLatticeSmoothDiracDelta3D(
+ const BlockLatticeSmoothDiracDelta3D& rhs)
+ :
+ BlockLatticeF3D(rhs._blockLattice, 3), _conv(rhs._conv), _cuboid(rhs._cuboid)
+{
+}
+
+template
+void BlockLatticeSmoothDiracDelta3D::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
--
cgit v1.2.3