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/boundary/extendedFiniteDifferenceBoundary2D.hh | 363 +++++++++++++++++++++
1 file changed, 363 insertions(+)
create mode 100644 src/boundary/extendedFiniteDifferenceBoundary2D.hh
(limited to 'src/boundary/extendedFiniteDifferenceBoundary2D.hh')
diff --git a/src/boundary/extendedFiniteDifferenceBoundary2D.hh b/src/boundary/extendedFiniteDifferenceBoundary2D.hh
new file mode 100644
index 0000000..70d1276
--- /dev/null
+++ b/src/boundary/extendedFiniteDifferenceBoundary2D.hh
@@ -0,0 +1,363 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2007 Orestis Malaspinas, Jonas Latt
+ * 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 EXTENDED_FINITE_DIFFERENCE_BOUNDARY_2D_HH
+#define EXTENDED_FINITE_DIFFERENCE_BOUNDARY_2D_HH
+
+#include "extendedFiniteDifferenceBoundary2D.h"
+#include "core/finiteDifference2D.h"
+#include "core/blockLattice2D.h"
+#include "core/util.h"
+#include "dynamics/lbHelpers.h"
+#include "dynamics/firstOrderLbHelpers.h"
+#include "boundaryInstantiator2D.h"
+
+
+namespace olb {
+
+
+/////////// ExtendedStraightFdBoundaryPostProcessor2D ///////////////////////////////////
+
+template
+ExtendedStraightFdBoundaryPostProcessor2D::
+ExtendedStraightFdBoundaryPostProcessor2D(int x0_, int x1_, int y0_, int y1_)
+ : x0(x0_), x1(x1_), y0(y0_), y1(y1_)
+{
+ OLB_PRECONDITION(x0==x1 || y0==y1);
+}
+
+template
+void ExtendedStraightFdBoundaryPostProcessor2D::
+processSubDomain(BlockLattice2D& blockLattice, int x0_, int x1_, int y0_, int y1_)
+{
+ using namespace olb::util::tensorIndices2D;
+ typedef lbHelpers lbH;
+ typedef DESCRIPTOR L;
+ enum {x,y};
+
+ int newX0, newX1, newY0, newY1;
+ if ( util::intersect (
+ x0, x1, y0, y1,
+ x0_, x1_, y0_, y1_,
+ newX0, newX1, newY0, newY1 ) ) {
+ for (int iX=newX0; iX<=newX1; ++iX) {
+ for (int iY=newY0; iY<=newY1; ++iY) {
+ Cell& cell = blockLattice.get(iX,iY);
+ T rho, u[L::d];
+ cell.computeRhoU(rho,u);
+
+ T uSqr = util::normSqr(u);
+
+ T dx_U[L::d], dy_U[L::d];
+ interpolateGradients<0>(blockLattice, dx_U, iX, iY);
+ interpolateGradients<1>(blockLattice, dy_U, iX, iY);
+
+ T rhoGradU[L::d][L::d];
+ rhoGradU[x][x] = rho * dx_U[x];
+ rhoGradU[x][y] = rho * dx_U[y];
+ rhoGradU[y][x] = rho * dy_U[x];
+ rhoGradU[y][y] = rho * dy_U[y];
+
+ T omega = blockLattice.getDynamics(iX, iY) -> getOmega();
+ T sToPi = - (T)1 / descriptors::invCs2() / omega;
+
+ T pi[util::TensorVal::n];
+ pi[xx] = (T)2 * rhoGradU[x][x] * sToPi;
+ pi[yy] = (T)2 * rhoGradU[y][y] * sToPi;
+ pi[xy] = (rhoGradU[x][y] + rhoGradU[y][x]) * sToPi;
+ // here ends the "regular" fdBoudaryCondition
+ // implemented in OpenLB
+
+ // first we compute the term
+ // (c_{i\alpha} \nabla_\beta)(rho*u_\alpha*u_\beta)
+ T dx_rho, dy_rho;
+ interpolateGradients<0>(blockLattice, dx_rho, iX, iY);
+ interpolateGradients<1>(blockLattice, dy_rho, iX, iY);
+ for (int iPop = 0; iPop < L::q; ++iPop) {
+ T cGradRhoUU = T();
+ for (int iAlpha=0; iAlpha < L::d; ++iAlpha) {
+ cGradRhoUU += descriptors::c(iPop,iAlpha) * (
+ dx_rho*u[iAlpha]*u[x] +
+ dx_U[iAlpha]*rho*u[x] +
+ dx_U[x]*rho*u[iAlpha] + //end of dx derivatice
+ dy_rho*u[iAlpha]*u[y] +
+ dy_U[iAlpha]*rho*u[y] +
+ dy_U[y]*rho*u[iAlpha]);
+ }
+
+ // then we compute the term
+ // c_{i\gamma}\nabla_{\gamma}(\rho*u_\alpha * u_\beta)
+ T cDivRhoUU[L::d][L::d]; //first step towards QcdivRhoUU
+ for (int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
+ for (int iBeta = 0; iBeta < L::d; ++iBeta) {
+ cDivRhoUU[iAlpha][iBeta] = descriptors::c(iPop,x) *
+ (dx_rho*u[iAlpha]*u[iBeta] +
+ dx_U[iAlpha]*rho*u[iBeta] +
+ dx_U[iBeta]*rho*u[iAlpha])
+ + descriptors::c(iPop,y) *
+ (dy_rho*u[iAlpha]*u[iBeta] +
+ dy_U[iAlpha]*rho*u[iBeta] +
+ dy_U[iBeta]*rho*u[iAlpha]);
+ }
+ }
+
+ //Finally we can compute
+ // Q_{i\alpha\beta}c_{i\gamma}\nabla_{\gamma}(\rho*u_\alpha * u_\beta)
+ // and Q_{i\alpha\beta}\rho\nabla_{\alpha}u_\beta
+ T qCdivRhoUU = T();
+ T qRhoGradU = T();
+ for (int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
+ for (int iBeta = 0; iBeta < L::d; ++iBeta) {
+ int ci_ci = descriptors::c(iPop,iAlpha)*descriptors::c(iPop,iBeta);
+ qCdivRhoUU += ci_ci * cDivRhoUU[iAlpha][iBeta];
+ qRhoGradU += ci_ci * rhoGradU[iAlpha][iBeta];
+ if (iAlpha == iBeta) {
+ qCdivRhoUU -= cDivRhoUU[iAlpha][iBeta]/descriptors::invCs2();
+ qRhoGradU -= rhoGradU[iAlpha][iBeta]/descriptors::invCs2();
+ }
+ }
+ }
+
+ // we then can reconstruct the value of the populations
+ // according to the complete C-E expansion term
+ cell[iPop] = lbH::equilibrium(iPop,rho,u,uSqr)
+ - descriptors::t(iPop) * descriptors::invCs2() / omega
+ * (qRhoGradU - cGradRhoUU + 0.5*descriptors::invCs2()*qCdivRhoUU);
+ }
+ }
+ }
+ }
+}
+
+template
+void ExtendedStraightFdBoundaryPostProcessor2D::
+process(BlockLattice2D& blockLattice)
+{
+ processSubDomain(blockLattice, x0, x1, y0, y1);
+}
+
+template
+template
+void ExtendedStraightFdBoundaryPostProcessor2D::
+interpolateGradients(BlockLattice2D const& blockLattice,
+ T velDeriv[DESCRIPTOR::d], int iX, int iY) const
+{
+ fd::DirectedGradients2D::
+ interpolateVector(velDeriv, blockLattice, iX, iY);
+}
+
+template
+template
+void ExtendedStraightFdBoundaryPostProcessor2D::
+interpolateGradients(BlockLattice2D const& blockLattice, T& rhoDeriv, int iX, int iY) const
+{
+ fd::DirectedGradients2D::
+ interpolateScalar(rhoDeriv, blockLattice, iX, iY);
+}
+
+
+//////// ExtendedStraightFdBoundaryPostProcessorGenerator ////////////////////////////////
+
+template
+ExtendedStraightFdBoundaryProcessorGenerator2D::
+ExtendedStraightFdBoundaryProcessorGenerator2D(int x0_, int x1_, int y0_, int y1_)
+ : PostProcessorGenerator2D(x0_, x1_, y0_, y1_)
+{ }
+
+template
+PostProcessor2D*
+ExtendedStraightFdBoundaryProcessorGenerator2D::generate() const
+{
+ return new ExtendedStraightFdBoundaryPostProcessor2D
+ ( this->x0, this->x1, this->y0, this->y1 );
+}
+
+template
+PostProcessorGenerator2D*
+ExtendedStraightFdBoundaryProcessorGenerator2D::clone() const
+{
+ return new ExtendedStraightFdBoundaryProcessorGenerator2D
+ (this->x0, this->x1, this->y0, this->y1);
+}
+
+
+////////// Class ExtendedFdBoundaryManager2D /////////////////////////////////////////
+
+template
+class ExtendedFdBoundaryManager2D {
+public:
+ template static Momenta*
+ getVelocityBoundaryMomenta();
+ template static Dynamics*
+ getVelocityBoundaryDynamics(T omega, Momenta& momenta);
+ template static PostProcessorGenerator2D*
+ getVelocityBoundaryProcessor(int x0, int x1, int y0, int y1);
+
+ template static Momenta*
+ getPressureBoundaryMomenta();
+ template static Dynamics*
+ getPressureBoundaryDynamics(T omega, Momenta& momenta);
+ template static PostProcessorGenerator2D*
+ getPressureBoundaryProcessor(int x0, int x1, int y0, int y1);
+
+ template static PostProcessorGenerator2D*
+ getConvectionBoundaryProcessor(int x0, int x1, int y0, int y1, T* uAv=NULL);
+
+ template static Momenta*
+ getExternalVelocityCornerMomenta();
+ template static Dynamics*
+ getExternalVelocityCornerDynamics(T omega, Momenta& momenta);
+ template static PostProcessorGenerator2D*
+ getExternalVelocityCornerProcessor(int x, int y);
+
+ template static Momenta*
+ getInternalVelocityCornerMomenta();
+ template static Dynamics*
+ getInternalVelocityCornerDynamics(T omega, Momenta& momenta);
+ template static PostProcessorGenerator2D*
+ getInternalVelocityCornerProcessor(int x, int y);
+};
+
+
+template
+template
+Momenta* ExtendedFdBoundaryManager2D::getVelocityBoundaryMomenta()
+{
+ return new BasicDirichletBM;
+}
+
+template
+template
+Dynamics* ExtendedFdBoundaryManager2D::
+getVelocityBoundaryDynamics(T omega, Momenta& momenta)
+{
+ return new MixinDynamics(omega, momenta);
+}
+
+template
+template
+PostProcessorGenerator2D* ExtendedFdBoundaryManager2D::
+getVelocityBoundaryProcessor(int x0, int x1, int y0, int y1)
+{
+ return new ExtendedStraightFdBoundaryProcessorGenerator2D
+ (x0, x1, y0, y1);
+}
+
+template
+template
+Momenta* ExtendedFdBoundaryManager2D::getPressureBoundaryMomenta()
+{
+ return new BasicDirichletBM;
+}
+
+template
+template
+Dynamics* ExtendedFdBoundaryManager2D::
+getPressureBoundaryDynamics(T omega, Momenta& momenta)
+{
+ return new MixinDynamics(omega, momenta);
+}
+
+template
+template
+PostProcessorGenerator2D*
+ExtendedFdBoundaryManager2D::
+getPressureBoundaryProcessor(int x0, int x1, int y0, int y1)
+{
+ return new ExtendedStraightFdBoundaryProcessorGenerator2D
+ (x0, x1, y0, y1);
+}
+
+template
+template
+PostProcessorGenerator2D*
+ExtendedFdBoundaryManager2D::
+getConvectionBoundaryProcessor(int x0, int x1, int y0, int y1, T* uAv)
+{
+ return nullptr;
+}
+
+template
+template
+Momenta*
+ExtendedFdBoundaryManager2D::getExternalVelocityCornerMomenta()
+{
+ return new FixedVelocityBM;
+}
+
+template
+template
+Dynamics* ExtendedFdBoundaryManager2D::
+getExternalVelocityCornerDynamics(T omega, Momenta& momenta)
+{
+ return new MixinDynamics(omega, momenta);
+}
+
+template
+template
+PostProcessorGenerator2D*
+ExtendedFdBoundaryManager2D::getExternalVelocityCornerProcessor(int x, int y)
+{
+ return new OuterVelocityCornerProcessorGenerator2D (x,y);
+}
+
+template
+template
+Momenta*
+ExtendedFdBoundaryManager2D::getInternalVelocityCornerMomenta()
+{
+ return new InnerCornerVelBM2D;
+}
+
+template
+template
+Dynamics* ExtendedFdBoundaryManager2D::
+getInternalVelocityCornerDynamics(T omega, Momenta& momenta)
+{
+ return new CombinedRLBdynamics(omega, momenta);
+}
+
+template
+template
+PostProcessorGenerator2D*
+ExtendedFdBoundaryManager2D::getInternalVelocityCornerProcessor (int x, int y)
+{
+ return nullptr;
+}
+
+
+////////// Factory functions //////////////////////////////////////////////////
+
+template
+OnLatticeBoundaryCondition2D*
+createExtendedFdBoundaryCondition2D(BlockLatticeStructure2D& block)
+{
+ return new BoundaryConditionInstantiator2D <
+ T, DESCRIPTOR,
+ ExtendedFdBoundaryManager2D > (block);
+}
+
+
+} // namespace olb
+
+#endif
--
cgit v1.2.3