From 8ba42a9a9dd2630a651d2cfa2e0ebc7fb1843100 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Wed, 1 May 2019 21:06:52 +0200 Subject: Move refinement to contrib --- global.mk | 2 + src/contrib/MakeHeader | 21 ++ src/contrib/contrib2D.h | 38 ++++ src/contrib/contrib2D.hh | 25 +++ src/contrib/module.mk | 27 +++ src/contrib/refinement/MakeHeader | 27 +++ src/contrib/refinement/coupler2D.cpp | 35 +++ src/contrib/refinement/coupler2D.h | 90 ++++++++ src/contrib/refinement/coupler2D.hh | 349 +++++++++++++++++++++++++++++ src/contrib/refinement/grid2D.cpp | 34 +++ src/contrib/refinement/grid2D.h | 169 ++++++++++++++ src/contrib/refinement/grid2D.hh | 388 +++++++++++++++++++++++++++++++++ src/contrib/refinement/module.mk | 27 +++ src/contrib/refinement/refinement2D.h | 27 +++ src/contrib/refinement/refinement2D.hh | 27 +++ src/olb2D.h | 1 - src/olb2D.hh | 1 - src/refinement/MakeHeader | 27 --- src/refinement/coupler2D.cpp | 35 --- src/refinement/coupler2D.h | 90 -------- src/refinement/coupler2D.hh | 339 ---------------------------- src/refinement/grid2D.cpp | 34 --- src/refinement/grid2D.h | 169 -------------- src/refinement/grid2D.hh | 388 --------------------------------- src/refinement/module.mk | 27 --- src/refinement/refinement2D.h | 27 --- src/refinement/refinement2D.hh | 27 --- 27 files changed, 1286 insertions(+), 1165 deletions(-) create mode 100644 src/contrib/MakeHeader create mode 100644 src/contrib/contrib2D.h create mode 100644 src/contrib/contrib2D.hh create mode 100644 src/contrib/module.mk create mode 100644 src/contrib/refinement/MakeHeader create mode 100644 src/contrib/refinement/coupler2D.cpp create mode 100644 src/contrib/refinement/coupler2D.h create mode 100644 src/contrib/refinement/coupler2D.hh create mode 100644 src/contrib/refinement/grid2D.cpp create mode 100644 src/contrib/refinement/grid2D.h create mode 100644 src/contrib/refinement/grid2D.hh create mode 100644 src/contrib/refinement/module.mk create mode 100644 src/contrib/refinement/refinement2D.h create mode 100644 src/contrib/refinement/refinement2D.hh delete mode 100644 src/refinement/MakeHeader delete mode 100644 src/refinement/coupler2D.cpp delete mode 100644 src/refinement/coupler2D.h delete mode 100644 src/refinement/coupler2D.hh delete mode 100644 src/refinement/grid2D.cpp delete mode 100644 src/refinement/grid2D.h delete mode 100644 src/refinement/grid2D.hh delete mode 100644 src/refinement/module.mk delete mode 100644 src/refinement/refinement2D.h delete mode 100644 src/refinement/refinement2D.hh diff --git a/global.mk b/global.mk index e2e8329..79cedc9 100644 --- a/global.mk +++ b/global.mk @@ -71,6 +71,8 @@ LIBS := -l$(LIB) -lz SUBDIRS := src/boundary \ src/communication \ src/dynamics \ + src/contrib \ + src/contrib/functors \ src/core \ src/geometry \ src/external/tinyxml \ diff --git a/src/contrib/MakeHeader b/src/contrib/MakeHeader new file mode 100644 index 0000000..67c5005 --- /dev/null +++ b/src/contrib/MakeHeader @@ -0,0 +1,21 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# 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. diff --git a/src/contrib/contrib2D.h b/src/contrib/contrib2D.h new file mode 100644 index 0000000..eeb5b79 --- /dev/null +++ b/src/contrib/contrib2D.h @@ -0,0 +1,38 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2007 the OpenLB project + * + * 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. +*/ + +/** \file + * Groups all the 2D include files in the contrib directory. + */ + +#include "advectionDiffusionVelocityCouplingPostProcessor2D.h" +#include "HerschelBulkleyBGKdynamics.h" +//#include "movingBoundaryBGKdynamics.h" +//#include "movingBoundaryFunctors.h" +#include "keepIncomingBoundary2D.h" +#include "keepIncomingDynamics.h" +//#include "straightBouzidiBoundary2D.h" +//#include "straightBouzidiDynamics.h" +#include "entropicLatticeDescriptors.h" +#include "optEntropicDynamics.h" +#include "shanLatticeDescriptors.h" +#include "SmagorinskyHerschelBulkleyPorousBGKdynamics.h" +//#include "graphLoadBalancer.h" +#include "refinement/refinement2D.h" diff --git a/src/contrib/contrib2D.hh b/src/contrib/contrib2D.hh new file mode 100644 index 0000000..00f6c20 --- /dev/null +++ b/src/contrib/contrib2D.hh @@ -0,0 +1,25 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2007 the OpenLB project + * + * 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. +*/ + +/** \file + * Groups all the generic 2D template files in the contrib directory. + */ + +#include "refinement/refinement2D.hh" diff --git a/src/contrib/module.mk b/src/contrib/module.mk new file mode 100644 index 0000000..6e42e3e --- /dev/null +++ b/src/contrib/module.mk @@ -0,0 +1,27 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2017 Markus Mohrhard +# 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. + +current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST))) + +include $(addsuffix /MakeHeader, $(current_dir)) + +LIB_OBJECTS += $(foreach file, $($(BUILDTYPE)), $(OBJDIR)/$(current_dir)$(file).o) 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 +# +# +# 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 + * + * + * 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>; +template class FineCoupler2D>; +template class CoarseCoupler2D>; + +} 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 + * + * + * 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 +class Coupler2D { +protected: + Grid2D& _coarse; + Grid2D& _fine; + + const int _coarseSize; + const int _fineSize; + const bool _vertical; + + const Vector _physOrigin; + + const Vector& getFineLatticeR(int y) const; + const Vector& getCoarseLatticeR(int y) const; + + T getScalingFactor() const; + T getInvScalingFactor() const; + +private: + std::vector> _coarseLatticeR; + std::vector> _fineLatticeR; + +public: + Coupler2D(Grid2D& coarse, Grid2D& fine, + Vector origin, Vector extend); + +}; + +template +class FineCoupler2D : public Coupler2D { +private: + std::vector> _c2f_rho; + std::vector> _c2f_u; + std::vector> _c2f_fneq; + +public: + FineCoupler2D(Grid2D& coarse, Grid2D& fine, + Vector origin, Vector extend); + + void store(); + void interpolate(); + void couple(); + +}; + +template +class CoarseCoupler2D : public Coupler2D { +public: + CoarseCoupler2D(Grid2D& coarse, Grid2D& fine, + Vector origin, Vector 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 + * + * + * 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 +T Coupler2D::getScalingFactor() const +{ + const T coarseTau = _coarse.getConverter().getLatticeRelaxationTime(); + return (coarseTau - 0.25) / coarseTau; +} + +template +T Coupler2D::getInvScalingFactor() const +{ + return 1./getScalingFactor(); +} + +template +const Vector& Coupler2D::getFineLatticeR(int y) const +{ + return _fineLatticeR[y]; +} + +template +const Vector& Coupler2D::getCoarseLatticeR(int y) const +{ + return _coarseLatticeR[y]; +} + +template +Coupler2D::Coupler2D(Grid2D& coarse, Grid2D& fine, + Vector origin, Vector 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 stepPhysR = _vertical ? Vector {0, deltaX} : + Vector {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 +FineCoupler2D::FineCoupler2D(Grid2D& coarse, Grid2D& fine, + Vector origin, Vector extend): + Coupler2D(coarse, fine, origin, extend), + _c2f_rho(this->_coarseSize), + _c2f_u(this->_coarseSize, Vector(T{})), + _c2f_fneq(this->_coarseSize, Vector(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 +void FineCoupler2D::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 coarseCell; + coarseLattice.get(pos, coarseCell); + lbHelpers::computeRhoU(coarseCell, rho, u); + lbHelpers::computeFneq(coarseCell, fNeq, rho, u); + + _c2f_rho[y] = Vector(rho); + _c2f_u[y] = Vector(u); + _c2f_fneq[y] = Vector(fNeq); + } +} + +template +Vector order2interpolation(const Vector& f0, const Vector& f1) +{ + return 0.5 * (f0 + f1); +} + +template +Vector order2interpolation(const std::vector>& data, int y) +{ + return 0.5 * (data[y] + data[y+1]); +} + +template +Vector order3interpolation(const std::vector>& 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 +Vector order4interpolation(const std::vector>& data, int y) +{ + return 9./16. * (data[y] + data[y+1]) - 1./16. * (data[y-1] + data[y+2]); +} + + +template +void FineCoupler2D::interpolate() +{ + auto& coarseLattice = this->_coarse.getSuperLattice(); + +#ifdef PARALLEL_MODE_OMP + #pragma omp parallel for +#endif + for (int y=0; y < this->_coarseSize; ++y) { + Cell coarseCell; + coarseLattice.get(this->getCoarseLatticeR(y), coarseCell); + + T rho{}; + T u[2] {}; + lbHelpers::computeRhoU(coarseCell, rho, u); + + _c2f_rho[y] = order2interpolation(Vector(rho), _c2f_rho[y]); + _c2f_u[y] = order2interpolation(Vector(u), _c2f_u[y]); + + T fNeq[DESCRIPTOR::q] {}; + lbHelpers::computeFneq(coarseCell, fNeq, rho, u); + + _c2f_fneq[y] = order2interpolation(Vector(fNeq), _c2f_fneq[y]); + } +} + +template +void FineCoupler2D::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 coarseCell; + coarseLattice.get(coarsePos, coarseCell); + lbHelpers::computeFeq(coarseCell, fEq); + + Cell 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 fineCell; + fineLattice.get(finePos, fineCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineCell[iPop] = lbHelpers::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 fineCell; + fineLattice.get(finePos, fineCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineCell[iPop] = lbHelpers::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 fineCell; + fineLattice.get(finePos, fineCell); + + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + fineCell[iPop] = lbHelpers::equilibrium(iPop, rho[0], u.data, uSqr) + + this->getScalingFactor() * fneq[iPop]; + } + + fineLattice.set(finePos, fineCell); + } +} + + +template +void computeRestrictedFneq(const SuperLattice2D& lattice, + Vector latticeR, + T restrictedFneq[DESCRIPTOR::q]) +{ + for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { + const auto neighbor = latticeR + Vector {0, descriptors::c(iPop,0), descriptors::c(iPop,1)}; + Cell cell; + lattice.get(neighbor, cell); + + T fNeq[DESCRIPTOR::q] {}; + lbHelpers::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 +CoarseCoupler2D::CoarseCoupler2D( + Grid2D& coarse, Grid2D& fine, + Vector origin, Vector extend): + Coupler2D(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 +void CoarseCoupler2D::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 fineCell; + fineLattice.get(finePos, fineCell); + lbHelpers::computeFeq(fineCell, fEq); + + T fNeq[DESCRIPTOR::q] {}; + computeRestrictedFneq(fineLattice, finePos, fNeq); + + Cell 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 + * + * + * 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>; +template class RefiningGrid2D>; + +} 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 + * + * + * 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 +#include +#include +#include + +#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 class FineCoupler2D; +template class CoarseCoupler2D; +template class RefiningGrid2D; + +template class sOnLatticeBoundaryCondition2D; +template class sOffLatticeBoundaryCondition2D; + +template +using RelaxationTime = utilities::NamedType; + +template +using LatticeVelocity = utilities::NamedType; + +template +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 +class Grid2D { +protected: + FunctorPtr> _domainF; + const Characteristics _characteristics; + + std::unique_ptr> _converter; + std::unique_ptr> _cuboids; + std::unique_ptr> _balancer; + std::unique_ptr> _geometry; + std::unique_ptr> _lattice; + + std::vector>> _dynamics; + std::vector>> _onLatticeBoundaryConditions; + std::vector>> _offLatticeBoundaryConditions; + + std::vector>> _fineGrids; + + std::vector>> _fineCouplers; + std::vector>> _coarseCouplers; + +public: + Grid2D(FunctorPtr>&& domainF, + RelaxationTime tau, + int resolution, + Characteristics characteristics); + Grid2D(FunctorPtr>&& domainF, + LatticeVelocity latticeVelocity, + int resolution, + Characteristics characteristics); + + Grid2D(FunctorPtr>&& domainF, RelaxationTime tau, int resolution, int re); + Grid2D(FunctorPtr>&& domainF, LatticeVelocity uMax, int resolution, int re); + + Characteristics getCharacteristics() const; + + UnitConverter& getConverter(); + CuboidGeometry2D& getCuboidGeometry(); + LoadBalancer& getLoadBalancer(); + SuperGeometry2D& getSuperGeometry(); + SuperLattice2D& getSuperLattice(); + + Dynamics& addDynamics(std::unique_ptr>&& dynamics); + sOnLatticeBoundaryCondition2D& getOnLatticeBoundaryCondition(); + sOffLatticeBoundaryCondition2D& getOffLatticeBoundaryCondition(); + + void collideAndStream(); + + FineCoupler2D& addFineCoupling( + Grid2D& fineGrid, Vector origin, Vector extend); + CoarseCoupler2D& addCoarseCoupling( + Grid2D& fineGrid, Vector origin, Vector extend); + + Vector alignOriginToGrid(Vector physR) const; + Vector alignExtendToGrid(Vector physR) const; + + RefiningGrid2D& refine(Vector origin, Vector extend, bool addCouplers=true); + + void forEachGrid(std::function&)>&& f); + void forEachGrid(const std::string& id, std::function&,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& locate(Vector pos); + + std::size_t getActiveVoxelN() const; + +}; + +template +class RefiningGrid2D : public Grid2D { +private: + const Vector _origin; + const Vector _extend; + + Grid2D& _parentGrid; + +public: + RefiningGrid2D(Grid2D& parentGrid, Vector origin, Vector extend); + + Vector getOrigin() const; + Vector getExtend() const; + + /// Indicates the subdomain of the coarse grid rendered moot by refinement + std::unique_ptr> 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 + * + * + * 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 +Grid2D::Grid2D(FunctorPtr>&& domainF, + RelaxationTime tau, + int resolution, + Characteristics characteristics): + _domainF(std::move(domainF)), + _characteristics(characteristics), + _converter(new UnitConverterFromResolutionAndRelaxationTime( + 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( + *_domainF, + _converter->getConversionFactorLength(), +#ifdef PARALLEL_MODE_MPI + singleton::mpi().getSize() +#else + 1 +#endif + )), + _balancer(new HeuristicLoadBalancer( + *_cuboids)), + _geometry(new SuperGeometry2D( + *_cuboids, + *_balancer, + 2)), + _lattice(new SuperLattice2D( + *_geometry)) +{ + _converter->print(); +} + +template +Grid2D::Grid2D(FunctorPtr>&& domainF, + LatticeVelocity latticeVelocity, + int resolution, + Characteristics characteristics): + _domainF(std::move(domainF)), + _characteristics(characteristics), + _converter(new UnitConverterFromResolutionAndLatticeVelocity( + 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( + *_domainF, + _converter->getConversionFactorLength(), +#ifdef PARALLEL_MODE_MPI + singleton::mpi().getSize() +#else + 1 +#endif + )), + _balancer(new HeuristicLoadBalancer( + *_cuboids)), + _geometry(new SuperGeometry2D( + *_cuboids, + *_balancer, + 2)), + _lattice(new SuperLattice2D( + *_geometry)) +{ + _converter->print(); +} + +template +Grid2D::Grid2D( + FunctorPtr>&& domainF, + RelaxationTime tau, + int resolution, + int re): + Grid2D(std::forward(domainF), + tau, + resolution, + Characteristics(re)) +{ } + +template +Grid2D::Grid2D( + FunctorPtr>&& domainF, + LatticeVelocity latticeVelocity, + int resolution, + int re): + Grid2D(std::forward(domainF), + latticeVelocity, + resolution, + Characteristics(re)) +{ } + +template +Characteristics Grid2D::getCharacteristics() const +{ + return _characteristics; +} + +template +UnitConverter& Grid2D::getConverter() +{ + return *_converter; +} + +template +CuboidGeometry2D& Grid2D::getCuboidGeometry() +{ + return *_cuboids; +} + +template +LoadBalancer& Grid2D::getLoadBalancer() +{ + return *_balancer; +} + +template +SuperGeometry2D& Grid2D::getSuperGeometry() +{ + return *_geometry; +} + +template +SuperLattice2D& Grid2D::getSuperLattice() +{ + return *_lattice; +} + +template +Dynamics& Grid2D::addDynamics( + std::unique_ptr>&& dynamics) +{ + Dynamics& ref = *dynamics; + _dynamics.emplace_back(std::move(dynamics)); + return ref; +} + +template +sOnLatticeBoundaryCondition2D& Grid2D::getOnLatticeBoundaryCondition() +{ + _onLatticeBoundaryConditions.emplace_back( + new sOnLatticeBoundaryCondition2D(getSuperLattice())); + return *_onLatticeBoundaryConditions.back(); +} + +template +sOffLatticeBoundaryCondition2D& Grid2D::getOffLatticeBoundaryCondition() +{ + _offLatticeBoundaryConditions.emplace_back( + new sOffLatticeBoundaryCondition2D(getSuperLattice())); + return *_offLatticeBoundaryConditions.back(); +} + +template +void Grid2D::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(); + } + + for ( auto& fineCoupler : _fineCouplers ) { + fineCoupler->store(); + fineCoupler->couple(); + } + + for ( auto& coarseCoupler : _coarseCouplers ) { + coarseCoupler->couple(); + } +} + +template +FineCoupler2D& Grid2D::addFineCoupling( + Grid2D& fineGrid, Vector origin, Vector extend) +{ + _fineCouplers.emplace_back( + new FineCoupler2D( + *this, fineGrid, origin, extend)); + return *_fineCouplers.back(); +} + +template +CoarseCoupler2D& Grid2D::addCoarseCoupling( + Grid2D& fineGrid, Vector origin, Vector extend) +{ + _coarseCouplers.emplace_back( + new CoarseCoupler2D( + *this, fineGrid, origin, extend)); + return *_coarseCouplers.back(); +} + +template +Vector Grid2D::alignOriginToGrid(Vector physR) const +{ + Vector latticeR{}; + _cuboids->getLatticeR(physR, latticeR); + return _cuboids->getPhysR(latticeR.toStdVector()); +} + +template +Vector Grid2D::alignExtendToGrid(Vector extend) const +{ + const T deltaX = _converter->getPhysDeltaX(); + return { + static_cast(std::floor(extend[0] / deltaX)) * deltaX, + static_cast(std::floor(extend[1] / deltaX)) * deltaX + }; +} + +template +void Grid2D::forEachGrid(std::function&)>&& f) +{ + f(*this); + for (auto& grid : _fineGrids) { + grid->forEachGrid(std::forward(f)); + } +} + +template +void Grid2D::forEachGrid( + const std::string& id, + std::function&,const std::string&)>&& f) +{ + f(*this, id); + for (std::size_t i = 0; i < _fineGrids.size(); ++i) { + _fineGrids[i]->forEachGrid(id + "_" + std::to_string(i), + std::forward(f)); + + } +} + +template +Grid2D& Grid2D::locate(Vector pos) +{ + int iC; + for (auto& grid : _fineGrids) { + if (grid->getCuboidGeometry().getC(pos, iC)) { + return grid->locate(pos); + } + } + return *this; +} + +template +std::size_t Grid2D::getActiveVoxelN() const +{ + std::size_t n = _geometry->getStatistics().getNvoxel(); + for (const auto& grid : _fineGrids) { + n += grid->getActiveVoxelN(); + } + return n; +} + +template +RefiningGrid2D& Grid2D::refine( + Vector wantedOrigin, Vector wantedExtend, bool addCouplers) +{ + if (addCouplers) { + auto& fineGrid = refine(wantedOrigin, wantedExtend, false); + + const Vector origin = fineGrid.getOrigin(); + const Vector extend = fineGrid.getExtend(); + + const Vector extendX = {extend[0],0}; + const Vector 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 innerOrigin = origin + coarseDeltaX; + const Vector innerExtendX = extendX - Vector {2*coarseDeltaX,0}; + const Vector innerExtendY = extendY - Vector {0,2*coarseDeltaX}; + + addCoarseCoupling(fineGrid, innerOrigin, innerExtendY); + addCoarseCoupling(fineGrid, innerOrigin + innerExtendX, innerExtendY); + addCoarseCoupling(fineGrid, innerOrigin + innerExtendY, innerExtendX); + addCoarseCoupling(fineGrid, innerOrigin, innerExtendX); + + return fineGrid; + } + else { + const Vector origin = alignOriginToGrid(wantedOrigin); + const Vector extend = alignExtendToGrid(wantedExtend); + + _fineGrids.emplace_back( + new RefiningGrid2D(*this, origin, extend)); + RefiningGrid2D& fineGrid = *_fineGrids.back(); + + auto refinedOverlap = fineGrid.getRefinedOverlap(); + _geometry->reset(*refinedOverlap); + + return fineGrid; + } +} + + +template +RefiningGrid2D::RefiningGrid2D( + Grid2D& parentGrid, Vector origin, Vector extend): + Grid2D( + std::unique_ptr>(new IndicatorCuboid2D(extend, origin)), + RelaxationTime(2*parentGrid.getConverter().getLatticeRelaxationTime() - 0.5), + 2*parentGrid.getConverter().getResolution(), + parentGrid.getCharacteristics()), + _origin(origin), + _extend(extend), + _parentGrid(parentGrid) { } + +template +Vector RefiningGrid2D::getOrigin() const +{ + return _origin; +} + +template +Vector RefiningGrid2D::getExtend() const +{ + return _extend; +} + +template +std::unique_ptr> RefiningGrid2D::getRefinedOverlap() const +{ + const T coarseDeltaX = _parentGrid.getConverter().getPhysDeltaX(); + + return std::unique_ptr>( + new IndicatorCuboid2D(_extend - 4*coarseDeltaX, _origin + 2*coarseDeltaX)); +} + +} + +#endif diff --git a/src/contrib/refinement/module.mk b/src/contrib/refinement/module.mk new file mode 100644 index 0000000..6e42e3e --- /dev/null +++ b/src/contrib/refinement/module.mk @@ -0,0 +1,27 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2017 Markus Mohrhard +# 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. + +current_dir := $(dir $(word $(words $(MAKEFILE_LIST)),$(MAKEFILE_LIST))) + +include $(addsuffix /MakeHeader, $(current_dir)) + +LIB_OBJECTS += $(foreach file, $($(BUILDTYPE)), $(OBJDIR)/$(current_dir)$(file).o) diff --git a/src/contrib/refinement/refinement2D.h b/src/contrib/refinement/refinement2D.h new file mode 100644 index 0000000..2dc726d --- /dev/null +++ b/src/contrib/refinement/refinement2D.h @@ -0,0 +1,27 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2019 Adrian Kummerländer + * + * 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_2D_H +#define REFINEMENT_2D_H + +#include "grid2D.h" +#include "coupler2D.h" + +#endif diff --git a/src/contrib/refinement/refinement2D.hh b/src/contrib/refinement/refinement2D.hh new file mode 100644 index 0000000..b21af48 --- /dev/null +++ b/src/contrib/refinement/refinement2D.hh @@ -0,0 +1,27 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2019 Adrian Kummerländer + * + * 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_2D_HH +#define REFINEMENT_2D_HH + +#include "grid2D.hh" +#include "coupler2D.hh" + +#endif diff --git a/src/olb2D.h b/src/olb2D.h index a5d72df..1f2caa1 100644 --- a/src/olb2D.h +++ b/src/olb2D.h @@ -7,4 +7,3 @@ #include #include #include -#include diff --git a/src/olb2D.hh b/src/olb2D.hh index 77189a9..481ab5b 100644 --- a/src/olb2D.hh +++ b/src/olb2D.hh @@ -7,4 +7,3 @@ #include #include #include -#include diff --git a/src/refinement/MakeHeader b/src/refinement/MakeHeader deleted file mode 100644 index bfcbf33..0000000 --- a/src/refinement/MakeHeader +++ /dev/null @@ -1,27 +0,0 @@ -# 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 -# -# -# 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/refinement/coupler2D.cpp b/src/refinement/coupler2D.cpp deleted file mode 100644 index 733d7c9..0000000 --- a/src/refinement/coupler2D.cpp +++ /dev/null @@ -1,35 +0,0 @@ -/* 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 - * - * - * 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>; -template class FineCoupler2D>; -template class CoarseCoupler2D>; - -} diff --git a/src/refinement/coupler2D.h b/src/refinement/coupler2D.h deleted file mode 100644 index 9b310b7..0000000 --- a/src/refinement/coupler2D.h +++ /dev/null @@ -1,90 +0,0 @@ -/* 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 - * - * - * 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 -class Coupler2D { -protected: - Grid2D& _coarse; - Grid2D& _fine; - - const int _coarseSize; - const int _fineSize; - const bool _vertical; - - const Vector _physOrigin; - - const Vector& getFineLatticeR(int y) const; - const Vector& getCoarseLatticeR(int y) const; - - T getScalingFactor() const; - T getInvScalingFactor() const; - -private: - std::vector> _coarseLatticeR; - std::vector> _fineLatticeR; - -public: - Coupler2D(Grid2D& coarse, Grid2D& fine, - Vector origin, Vector extend); - -}; - -template -class FineCoupler2D : public Coupler2D { -private: - std::vector> _c2f_rho; - std::vector> _c2f_u; - std::vector> _c2f_fneq; - -public: - FineCoupler2D(Grid2D& coarse, Grid2D& fine, - Vector origin, Vector extend); - - void store(); - void interpolate(); - void couple(); - -}; - -template -class CoarseCoupler2D : public Coupler2D { -public: - CoarseCoupler2D(Grid2D& coarse, Grid2D& fine, - Vector origin, Vector extend); - - void couple(); - -}; - - -} - -#endif diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh deleted file mode 100644 index 800986d..0000000 --- a/src/refinement/coupler2D.hh +++ /dev/null @@ -1,339 +0,0 @@ -/* 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 - * - * - * 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 -T Coupler2D::getScalingFactor() const -{ - const T coarseTau = _coarse.getConverter().getLatticeRelaxationTime(); - return (coarseTau - 0.25) / coarseTau; -} - -template -T Coupler2D::getInvScalingFactor() const -{ - return 1./getScalingFactor(); -} - -template -const Vector& Coupler2D::getFineLatticeR(int y) const -{ - return _fineLatticeR[y]; -} - -template -const Vector& Coupler2D::getCoarseLatticeR(int y) const -{ - return _coarseLatticeR[y]; -} - -template -Coupler2D::Coupler2D(Grid2D& coarse, Grid2D& fine, - Vector origin, Vector 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 stepPhysR = _vertical ? Vector {0, deltaX} : - Vector {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 -FineCoupler2D::FineCoupler2D(Grid2D& coarse, Grid2D& fine, - Vector origin, Vector extend): - Coupler2D(coarse, fine, origin, extend), - _c2f_rho(this->_coarseSize), - _c2f_u(this->_coarseSize, Vector(T{})), - _c2f_fneq(this->_coarseSize, Vector(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 -void FineCoupler2D::store() -{ - auto& coarseLattice = this->_coarse.getSuperLattice(); - - #pragma omp parallel for - for (int y=0; y < this->_coarseSize; ++y) { - const auto pos = this->getCoarseLatticeR(y); - T rho{}; - T u[2] {}; - T fNeq[DESCRIPTOR::q] {}; - Cell coarseCell; - coarseLattice.get(pos, coarseCell); - lbHelpers::computeRhoU(coarseCell, rho, u); - lbHelpers::computeFneq(coarseCell, fNeq, rho, u); - - _c2f_rho[y] = Vector(rho); - _c2f_u[y] = Vector(u); - _c2f_fneq[y] = Vector(fNeq); - } -} - -template -Vector order2interpolation(const Vector& f0, const Vector& f1) -{ - return 0.5 * (f0 + f1); -} - -template -Vector order2interpolation(const std::vector>& data, int y) -{ - return 0.5 * (data[y] + data[y+1]); -} - -template -Vector order3interpolation(const std::vector>& 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 -Vector order4interpolation(const std::vector>& data, int y) -{ - return 9./16. * (data[y] + data[y+1]) - 1./16. * (data[y-1] + data[y+2]); -} - - -template -void FineCoupler2D::interpolate() -{ - auto& coarseLattice = this->_coarse.getSuperLattice(); - - #pragma omp parallel for - for (int y=0; y < this->_coarseSize; ++y) { - Cell coarseCell; - coarseLattice.get(this->getCoarseLatticeR(y), coarseCell); - - T rho{}; - T u[2] {}; - lbHelpers::computeRhoU(coarseCell, rho, u); - - _c2f_rho[y] = order2interpolation(Vector(rho), _c2f_rho[y]); - _c2f_u[y] = order2interpolation(Vector(u), _c2f_u[y]); - - T fNeq[DESCRIPTOR::q] {}; - lbHelpers::computeFneq(coarseCell, fNeq, rho, u); - - _c2f_fneq[y] = order2interpolation(Vector(fNeq), _c2f_fneq[y]); - } -} - -template -void FineCoupler2D::couple() -{ - const auto& coarseLattice = this->_coarse.getSuperLattice(); - auto& fineLattice = this->_fine.getSuperLattice(); - - #pragma omp parallel for - 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 coarseCell; - coarseLattice.get(coarsePos, coarseCell); - lbHelpers::computeFeq(coarseCell, fEq); - - Cell 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); - } - - #pragma omp parallel for - 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 fineCell; - fineLattice.get(finePos, fineCell); - - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - fineCell[iPop] = lbHelpers::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 fineCell; - fineLattice.get(finePos, fineCell); - - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - fineCell[iPop] = lbHelpers::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 fineCell; - fineLattice.get(finePos, fineCell); - - for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) { - fineCell[iPop] = lbHelpers::equilibrium(iPop, rho[0], u.data, uSqr) -