diff options
Diffstat (limited to 'src/contrib/refinement')
| -rw-r--r-- | src/contrib/refinement/MakeHeader | 27 | ||||
| -rw-r--r-- | src/contrib/refinement/coupler2D.cpp | 35 | ||||
| -rw-r--r-- | src/contrib/refinement/coupler2D.h | 90 | ||||
| -rw-r--r-- | src/contrib/refinement/coupler2D.hh | 349 | ||||
| -rw-r--r-- | src/contrib/refinement/grid2D.cpp | 34 | ||||
| -rw-r--r-- | src/contrib/refinement/grid2D.h | 169 | ||||
| -rw-r--r-- | src/contrib/refinement/grid2D.hh | 388 | ||||
| -rw-r--r-- | src/contrib/refinement/module.mk | 27 | ||||
| -rw-r--r-- | src/contrib/refinement/refinement2D.h | 27 | ||||
| -rw-r--r-- | src/contrib/refinement/refinement2D.hh | 27 |
10 files changed, 1173 insertions, 0 deletions
diff --git a/src/contrib/refinement/MakeHeader b/src/contrib/refinement/MakeHeader new file mode 100644 index 0000000..bfcbf33 --- /dev/null +++ b/src/contrib/refinement/MakeHeader @@ -0,0 +1,27 @@ +# 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. + + +generic := + +precompiled := grid2D \ + coupler2D diff --git a/src/contrib/refinement/coupler2D.cpp b/src/contrib/refinement/coupler2D.cpp new file mode 100644 index 0000000..733d7c9 --- /dev/null +++ b/src/contrib/refinement/coupler2D.cpp @@ -0,0 +1,35 @@ +/* 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. + */ + +#include "coupler2D.h" +#include "coupler2D.hh" + +#include "dynamics/latticeDescriptors.h" + +namespace olb { + +template class Coupler2D<double,descriptors::D2Q9<>>; +template class FineCoupler2D<double,descriptors::D2Q9<>>; +template class CoarseCoupler2D<double,descriptors::D2Q9<>>; + +} diff --git a/src/contrib/refinement/coupler2D.h b/src/contrib/refinement/coupler2D.h new file mode 100644 index 0000000..9b310b7 --- /dev/null +++ b/src/contrib/refinement/coupler2D.h @@ -0,0 +1,90 @@ +/* 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_COUPLER_2D_H +#define REFINEMENT_COUPLER_2D_H + +#include "grid2D.h" + +namespace olb { + + +template <typename T, typename DESCRIPTOR> +class Coupler2D { +protected: + Grid2D<T,DESCRIPTOR>& _coarse; + Grid2D<T,DESCRIPTOR>& _fine; + + const int _coarseSize; + const int _fineSize; + const bool _vertical; + + const Vector<T,2> _physOrigin; + + const Vector<int,3>& getFineLatticeR(int y) const; + const Vector<int,3>& getCoarseLatticeR(int y) const; + + T getScalingFactor() const; + T getInvScalingFactor() const; + +private: + std::vector<Vector<int,3>> _coarseLatticeR; + std::vector<Vector<int,3>> _fineLatticeR; + +public: + Coupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, + Vector<T,2> origin, Vector<T,2> extend); + +}; + +template <typename T, typename DESCRIPTOR> +class FineCoupler2D : public Coupler2D<T,DESCRIPTOR> { +private: + std::vector<Vector<T,1>> _c2f_rho; + std::vector<Vector<T,DESCRIPTOR::d>> _c2f_u; + std::vector<Vector<T,DESCRIPTOR::q>> _c2f_fneq; + +public: + FineCoupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, + Vector<T,2> origin, Vector<T,2> extend); + + void store(); + void interpolate(); + void couple(); + +}; + +template <typename T, typename DESCRIPTOR> +class CoarseCoupler2D : public Coupler2D<T,DESCRIPTOR> { +public: + CoarseCoupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, + Vector<T,2> origin, Vector<T,2> extend); + + void couple(); + +}; + + +} + +#endif diff --git a/src/contrib/refinement/coupler2D.hh b/src/contrib/refinement/coupler2D.hh new file mode 100644 index 0000000..4c8b424 --- /dev/null +++ b/src/contrib/refinement/coupler2D.hh @@ -0,0 +1,349 @@ +/* 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_COUPLER_2D_HH +#define REFINEMENT_COUPLER_2D_HH + +#include "coupler2D.h" + +#include "dynamics/lbHelpers.h" + +namespace olb { + +template <typename T, typename DESCRIPTOR> +T Coupler2D<T,DESCRIPTOR>::getScalingFactor() const +{ + const T coarseTau = _coarse.getConverter().getLatticeRelaxationTime(); + return (coarseTau - 0.25) / coarseTau; +} + +template <typename T, typename DESCRIPTOR> +T Coupler2D<T,DESCRIPTOR>::getInvScalingFactor() const +{ + return 1./getScalingFactor(); +} + +template <typename T, typename DESCRIPTOR> +const Vector<int,3>& Coupler2D<T,DESCRIPTOR>::getFineLatticeR(int y) const +{ + return _fineLatticeR[y]; +} + +template <typename T, typename DESCRIPTOR> +const Vector<int,3>& Coupler2D<T,DESCRIPTOR>::getCoarseLatticeR(int y) const +{ + return _coarseLatticeR[y]; +} + +template <typename T, typename DESCRIPTOR> +Coupler2D<T,DESCRIPTOR>::Coupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, + Vector<T,2> origin, Vector<T,2> extend): + _coarse(coarse), + _fine(fine), + _coarseSize(floor(extend.norm() / coarse.getConverter().getPhysDeltaX() + 0.5)+1), + _fineSize(2*_coarseSize-1), + _vertical(util::nearZero(extend[0])), + _physOrigin(_coarse.alignOriginToGrid(origin)), + _coarseLatticeR(_coarseSize), + _fineLatticeR(_fineSize) +{ + OLB_ASSERT(util::nearZero(extend[0]) || util::nearZero(extend[1]), "Coupling is only implemented alongside unit vectors"); + + const auto& coarseGeometry = _coarse.getCuboidGeometry(); + const auto& fineGeometry = _fine.getCuboidGeometry(); + + const T deltaX = _fine.getConverter().getPhysDeltaX(); + const Vector<T,2> stepPhysR = _vertical ? Vector<T,2> {0, deltaX} : + Vector<T,2> {deltaX, 0}; + + for (int i=0; i < _fineSize; ++i) { + if (i % 2 == 0) { + coarseGeometry.getLatticeR(_physOrigin + i*stepPhysR, _coarseLatticeR[i/2]); + } + + fineGeometry.getLatticeR(_physOrigin + i*stepPhysR, _fineLatticeR[i]); + } +} + + +template <typename T, typename DESCRIPTOR> +FineCoupler2D<T,DESCRIPTOR>::FineCoupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, + Vector<T,2> origin, Vector<T,2> extend): + Coupler2D<T,DESCRIPTOR>(coarse, fine, origin, extend), + _c2f_rho(this->_coarseSize), + _c2f_u(this->_coarseSize, Vector<T,2>(T{})), + _c2f_fneq(this->_coarseSize, Vector<T,DESCRIPTOR::q>(T{})) +{ + OstreamManager clout(std::cout,"C2F"); + + const auto& coarseOrigin = this->getCoarseLatticeR(0); + const auto& fineOrigin = this->getFineLatticeR(0); + + clout << "coarse origin: " << coarseOrigin[0] << " " << coarseOrigin[1] << " " << coarseOrigin[2] << std::endl; + clout << "fine origin: " << fineOrigin[0] << " " << fineOrigin[1] << " " << fineOrigin[2] << std::endl; + clout << "fine size: " << this->_fineSize << std::endl; +} + +template <typename T, typename DESCRIPTOR> +void FineCoupler2D<T,DESCRIPTOR>::store() +{ + auto& coarseLattice = this->_coarse.getSuperLattice(); + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (int y=0; y < this->_coarseSize; ++y) { + const auto pos = this->getCoarseLatticeR(y); + T rho{}; + T u[2] {}; + T fNeq[DESCRIPTOR::q] {}; + Cell<T,DESCRIPTOR> coarseCell; + coarseLattice.get(pos, coarseCell); + lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, rho, u); + lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u); + + _c2f_rho[y] = Vector<T,1>(rho); + _c2f_u[y] = Vector<T,2>(u); + _c2f_fneq[y] = Vector<T,DESCRIPTOR::q>(fNeq); + } +} + +template <typename T, unsigned N> +Vector<T,N> order2interpolation(const Vector<T,N>& f0, const Vector<T,N>& f1) +{ + return 0.5 * (f0 + f1); +} + +template <typename T, unsigned N> +Vector<T,N> order2interpolation(const std::vector<Vector<T,N>>& data, int y) +{ + return 0.5 * (data[y] + data[y+1]); +} + +template <typename T, unsigned N> +Vector<T,N> order3interpolation(const std::vector<Vector<T,N>>& data, int y, bool ascending) +{ + if (ascending) { + return 3./8. * data[y] + 3./4. * data[y+1] - 1./8. * data[y+2]; + } + else { + return 3./8. * data[y] + 3./4. * data[y-1] - 1./8. * data[y-2]; + } +} + +template <typename T, unsigned N> +Vector<T,N> order4interpolation(const std::vector<Vector<T,N>>& data, int y) +{ + return 9./16. * (data[y] + data[y+1]) - 1./16. * (data[y-1] + data[y+2]); +} + + +template <typename T, typename DESCRIPTOR> +void FineCoupler2D<T,DESCRIPTOR>::interpolate() +{ + auto& coarseLattice = this->_coarse.getSuperLattice(); + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (int y=0; y < this->_coarseSize; ++y) { + Cell<T,DESCRIPTOR> coarseCell; + coarseLattice.get(this->getCoarseLatticeR(y), coarseCell); + + T rho{}; + T u[2] {}; + lbHelpers<T,DESCRIPTOR>::computeRhoU(coarseCell, rho, u); + + _c2f_rho[y] = order2interpolation(Vector<T,1>(rho), _c2f_rho[y]); + _c2f_u[y] = order2interpolation(Vector<T,2>(u), _c2f_u[y]); + + T fNeq[DESCRIPTOR::q] {}; + lbHelpers<T,DESCRIPTOR>::computeFneq(coarseCell, fNeq, rho, u); + + _c2f_fneq[y] = order2interpolation(Vector<T,DESCRIPTOR::q>(fNeq), _c2f_fneq[y]); + } +} + +template <typename T, typename DESCRIPTOR> +void FineCoupler2D<T,DESCRIPTOR>::couple() +{ + const auto& coarseLattice = this->_coarse.getSuperLattice(); + auto& fineLattice = this->_fine.getSuperLattice(); + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (int y=0; y < this->_coarseSize; ++y) { + const auto& coarsePos = this->getCoarseLatticeR(y); + const auto& finePos = this->getFineLatticeR(2*y); + + T fEq[DESCRIPTOR::q] {}; + Cell<T,DESCRIPTOR> coarseCell; + coarseLattice.get(coarsePos, coarseCell); + lbHelpers<T,DESCRIPTOR>::computeFeq(coarseCell, fEq); + + Cell<T,DESCRIPTOR> cell; + fineLattice.get(finePos, cell); + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + cell[iPop] = fEq[iPop] + this->getScalingFactor() * _c2f_fneq[y][iPop]; + } + fineLattice.set(finePos, cell); + } + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (int y=1; y < this->_coarseSize-2; ++y) { + const auto rho = order4interpolation(_c2f_rho, y); + const auto u = order4interpolation(_c2f_u, y); + const auto fneq = order4interpolation(_c2f_fneq, y); + + const T uSqr = u*u; + + const auto finePos = this->getFineLatticeR(1+2*y); + Cell<T,DESCRIPTOR> fineCell; + fineLattice.get(finePos, fineCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + + this->getScalingFactor() * fneq[iPop]; + } + + fineLattice.set(finePos, fineCell); + } + + { + const auto rho = order3interpolation(_c2f_rho, 0, true); + const auto u = order3interpolation(_c2f_u, 0, true); + const auto fneq = order3interpolation(_c2f_fneq, 0, true); + + const T uSqr = u*u; + + const auto& finePos = this->getFineLatticeR(1); + Cell<T,DESCRIPTOR> fineCell; + fineLattice.get(finePos, fineCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + + this->getScalingFactor() * fneq[iPop]; + } + + fineLattice.set(finePos, fineCell); + } + + { + const auto rho = order3interpolation(_c2f_rho, this->_coarseSize-1, false); + const auto u = order3interpolation(_c2f_u, this->_coarseSize-1, false); + const auto fneq = order3interpolation(_c2f_fneq, this->_coarseSize-1, false); + + const T uSqr = u*u; + + const auto& finePos = this->getFineLatticeR(this->_fineSize-2); + Cell<T,DESCRIPTOR> fineCell; + fineLattice.get(finePos, fineCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineCell[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho[0], u.data, uSqr) + + this->getScalingFactor() * fneq[iPop]; + } + + fineLattice.set(finePos, fineCell); + } +} + + +template <typename T, typename DESCRIPTOR> +void computeRestrictedFneq(const SuperLattice2D<T,DESCRIPTOR>& lattice, + Vector<int,3> latticeR, + T restrictedFneq[DESCRIPTOR::q]) +{ + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + const auto neighbor = latticeR + Vector<int,3> {0, descriptors::c<DESCRIPTOR>(iPop,0), descriptors::c<DESCRIPTOR>(iPop,1)}; + Cell<T,DESCRIPTOR> cell; + lattice.get(neighbor, cell); + + T fNeq[DESCRIPTOR::q] {}; + lbHelpers<T,DESCRIPTOR>::computeFneq(cell, fNeq); + + for (int jPop=0; jPop < DESCRIPTOR::q; ++jPop) { + restrictedFneq[jPop] += fNeq[jPop]; + } + } + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + restrictedFneq[iPop] /= DESCRIPTOR::q; + } +} + +template <typename T, typename DESCRIPTOR> +CoarseCoupler2D<T,DESCRIPTOR>::CoarseCoupler2D( + Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine, + Vector<T,2> origin, Vector<T,2> extend): + Coupler2D<T,DESCRIPTOR>(coarse, fine, origin, extend) +{ + OstreamManager clout(std::cout,"F2C"); + + const auto& coarseOrigin = this->getCoarseLatticeR(0); + const auto& fineOrigin = this->getFineLatticeR(0); + + clout << "coarse origin: " << coarseOrigin[0] << " " << coarseOrigin[1] << " " << coarseOrigin[2] << std::endl; + clout << "fine origin: " << fineOrigin[0] << " " << fineOrigin[1] << " " << fineOrigin[2] << std::endl; + clout << "coarse size: " << this->_coarseSize << std::endl; +} + +template <typename T, typename DESCRIPTOR> +void CoarseCoupler2D<T,DESCRIPTOR>::couple() +{ + const auto& fineLattice = this->_fine.getSuperLattice(); + auto& coarseLattice = this->_coarse.getSuperLattice(); + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (int y=0; y < this->_coarseSize; ++y) { + const auto& finePos = this->getFineLatticeR(2*y); + const auto& coarsePos = this->getCoarseLatticeR(y); + + T fEq[DESCRIPTOR::q] {}; + Cell<T,DESCRIPTOR> fineCell; + fineLattice.get(finePos, fineCell); + lbHelpers<T,DESCRIPTOR>::computeFeq(fineCell, fEq); + + T fNeq[DESCRIPTOR::q] {}; + computeRestrictedFneq(fineLattice, finePos, fNeq); + + Cell<T,DESCRIPTOR> coarseCell; + coarseLattice.get(coarsePos, coarseCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + coarseCell[iPop] = fEq[iPop] + this->getInvScalingFactor() * fNeq[iPop]; + } + + coarseLattice.set(coarsePos, coarseCell); + } +} + + +} + +#endif diff --git a/src/contrib/refinement/grid2D.cpp b/src/contrib/refinement/grid2D.cpp new file mode 100644 index 0000000..4472fc7 --- /dev/null +++ b/src/contrib/refinement/grid2D.cpp @@ -0,0 +1,34 @@ +/* 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. + */ + +#include "grid2D.h" +#include "grid2D.hh" + +#include "dynamics/latticeDescriptors.h" + +namespace olb { + +template class Grid2D<double,descriptors::D2Q9<>>; +template class RefiningGrid2D<double,descriptors::D2Q9<>>; + +} diff --git a/src/contrib/refinement/grid2D.h b/src/contrib/refinement/grid2D.h new file mode 100644 index 0000000..51d19ce --- /dev/null +++ b/src/contrib/refinement/grid2D.h @@ -0,0 +1,169 @@ +/* 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_H +#define REFINEMENT_GRID_2D_H + +#include <memory> +#include <vector> +#include <string> +#include <functional> + +#include "core/unitConverter.h" +#include "core/superLattice2D.h" +#include "geometry/superGeometry2D.h" +#include "geometry/cuboidGeometry2D.h" +#include "geometry/superGeometry2D.h" +#include "communication/loadBalancer.h" +#include "functors/analytical/indicator/indicatorF2D.h" +#include "utilities/functorPtr.h" +#include "utilities/namedType.h" + +namespace olb { + + +template <typename T, typename DESCRIPTOR> class FineCoupler2D; +template <typename T, typename DESCRIPTOR> class CoarseCoupler2D; +template <typename T, typename DESCRIPTOR> class RefiningGrid2D; + +template <typename T, typename DESCRIPTOR> class sOnLatticeBoundaryCondition2D; +template <typename T, typename DESCRIPTOR> class sOffLatticeBoundaryCondition2D; + +template <typename T> +using RelaxationTime = utilities::NamedType<T,struct NamedRelaxationTime>; + +template <typename T> +using LatticeVelocity = utilities::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, typename 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; + std::unique_ptr<LoadBalancer<T>> _balancer; + std::unique_ptr<SuperGeometry2D<T>> _geometry; + std::unique_ptr<SuperLattice2D<T,DESCRIPTOR>> _lattice; + + 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; + + std::vector<std::unique_ptr<FineCoupler2D<T,DESCRIPTOR>>> _fineCouplers; + std::vector<std::unique_ptr<CoarseCoupler2D<T,DESCRIPTOR>>> _coarseCouplers; + +public: + 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(); + LoadBalancer<T>& getLoadBalancer(); + SuperGeometry2D<T>& getSuperGeometry(); + SuperLattice2D<T,DESCRIPTOR>& getSuperLattice(); + + Dynamics<T,DESCRIPTOR>& addDynamics(std::unique_ptr<Dynamics<T,DESCRIPTOR>>&& dynamics); + sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& getOnLatticeBoundaryCondition(); + sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& getOffLatticeBoundaryCondition(); + + void collideAndStream(); + + FineCoupler2D<T,DESCRIPTOR>& addFineCoupling( + Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend); + CoarseCoupler2D<T,DESCRIPTOR>& addCoarseCoupling( + Grid2D<T,DESCRIPTOR>& fineGrid, Vector<T,2> origin, Vector<T,2> extend); + + Vector<T,2> alignOriginToGrid(Vector<T,2> physR) const; + Vector<T,2> alignExtendToGrid(Vector<T,2> physR) const; + + RefiningGrid2D<T,DESCRIPTOR>& refine(Vector<T,2> origin, Vector<T,2> extend, bool addCouplers=true); + + void forEachGrid(std::function<void(Grid2D<T,DESCRIPTOR>&)>&& f); + void forEachGrid(const std::string& id, std::function<void(Grid2D<T,DESCRIPTOR>&,const std::string&)>&& f); + + /// Returns the finest grid representing a physical position + /** + * Only works if pos is actually contained in a node of the refinement tree. + **/ + Grid2D<T,DESCRIPTOR>& locate(Vector<T,2> pos); + + std::size_t getActiveVoxelN() const; + +}; + +template <typename T, typename DESCRIPTOR> +class RefiningGrid2D : public Grid2D<T,DESCRIPTOR> { +private: + const Vector<T,2> _origin; + const Vector<T,2> _extend; + + Grid2D<T,DESCRIPTOR>& _parentGrid; + +public: + RefiningGrid2D(Grid2D<T,DESCRIPTOR>& parentGrid, Vector<T,2> origin, Vector<T,2> extend); + + Vector<T,2> getOrigin() const; + Vector<T,2> getExtend() const; + + /// Indicates the subdomain of the coarse grid rendered moot by refinement + std::unique_ptr<IndicatorF2D<T>> getRefinedOverlap() const; + +}; + + +} + +#endif diff --git a/src/contrib/refinement/grid2D.hh b/src/contrib/refinement/grid2D.hh new file mode 100644 index 0000000..d4310e5 --- /dev/null +++ b/src/contrib/refinement/grid2D.hh @@ -0,0 +1,388 @@ +/* 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 "utilities/vectorHelpers.h" +#include "boundary/superBoundaryCondition2D.h" +#include "boundary/superOffBoundaryCondition2D.h" +#include "communication/heuristicLoadBalancer.h" + +namespace olb { + + +template <typename T, typename DESCRIPTOR> +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! + 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(), +#ifdef PARALLEL_MODE_MPI + singleton::mpi().getSize() +#else + 1 +#endif + )), + _balancer(new HeuristicLoadBalancer<T>( + *_cuboids)), + _geometry(new SuperGeometry2D<T>( + *_cuboids, + *_balancer, + 2)), + _lattice(new SuperLattice2D<T,DESCRIPTOR>( + *_geometry)) +{ + _converter->print(); +} + +template <typename T, typename DESCRIPTOR> +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 + 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(), +#ifdef PARALLEL_MODE_MPI + singleton::mpi().getSize() +#else + 1 +#endif + )), + _balancer(new HeuristicLoadBalancer<T>( + *_cuboids)), + _geometry(new SuperGeometry2D<T>( + *_cuboids, + *_balancer, + 2)), + _lattice(new SuperLattice2D<T,DESCRIPTOR>( + *_geometry)) +{ + _converter->print(); +} + +template <typename T, typename 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, typename 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, typename DESCRIPTOR> +Characteristics<T> Grid2D<T,DESCRIPTOR>::getCharacteristics() const +{ + return _characteristics; +} + +template <typename T, typename DESCRIPTOR> +UnitConverter<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getConverter() +{ + return *_converter; +} + +template <typename T, typename DESCRIPTOR> +CuboidGeometry2D<T>& Grid2D<T,DESCRIPTOR>::getCuboidGeometry() +{ + return *_cuboids; +} + +template <typename T, typename DESCRIPTOR> +LoadBalancer<T>& Grid2D<T,DESCRIPTOR>::getLoadBalancer() +{ + return *_balancer; +} + +template <typename T, typename DESCRIPTOR> +SuperGeometry2D<T>& Grid2D<T,DESCRIPTOR>::getSuperGeometry() +{ + return *_geometry; +} + +template <typename T, typename DESCRIPTOR> +SuperLattice2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getSuperLattice() +{ + return *_lattice; +} + +template <typename T, typename DESCRIPTOR> +Dynamics<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::addDynamics( + std::unique_ptr<Dynamics<T,DESCRIPTOR>>&& dynamics) +{ + Dynamics<T,DESCRIPTOR>& ref = *dynamics; + _dynamics.emplace_back(std::move(dynamics)); + return ref; +} + +template <typename T, typename DESCRIPTOR> +sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getOnLatticeBoundaryCondition() +{ + _onLatticeBoundaryConditions.emplace_back( + new sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>(getSuperLattice())); + return *_onLatticeBoundaryConditions.back(); +} + +template <typename T, typename DESCRIPTOR> +sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& Grid2D<T,DESCRIPTOR>::getOffLatticeBoundaryCondition() +{ + _offLatticeBoundaryConditions.emplace_back( + new sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>(getSuperLattice())); + return *_offLatticeBoundaryConditions.back(); +} + +template <typename T, typename DESCRIPTOR> +void Grid2D<T,DESCRIPTOR>::collideAndStream() +{ + for ( auto& fineCoupler : _fineCouplers ) { + fineCoupler->store(); + } + + this->getSuperLattice().collideAndStream(); + + for ( auto& fineGrid : _fineGrids ) { + fineGrid->collideAndStream(); + } + + for ( auto& fineCoupler : _fineCouplers ) { + fineCoupler->interpolate(); + fineCoupler->couple(); + } + + for ( auto& fineGrid : _fineGrids ) { + fineGrid->collideAndStream( |