1 files changed, 136 insertions, 0 deletions
diff --git a/src/boundary/zouHeDynamics.hh b/src/boundary/zouHeDynamics.hh
new file mode 100644
index 0000000..cf3af15
--- /dev/null
+++ b/src/boundary/zouHeDynamics.hh
@@ -0,0 +1,136 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2006,2007 Orestis Malaspinas and Jonas Latt
+ *  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 ZOU_HE_DYNAMICS_HH
+#define ZOU_HE_DYNAMICS_HH
+
+#include "zouHeDynamics.h"
+#include "dynamics/latticeDescriptors.h"
+#include "core/util.h"
+#include "dynamics/lbHelpers.h"
+#include <cmath>
+
+
+namespace olb {
+
+ 
+
+template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation>
+ZouHeDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>::ZouHeDynamics (
+  T omega_, Momenta<T,DESCRIPTOR>& momenta_ )
+  : BasicDynamics<T,DESCRIPTOR>(momenta_),
+    boundaryDynamics(omega_, momenta_)
+{ }
+
+template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation>
+T ZouHeDynamics<T,DESCRIPTOR, Dynamics, direction, orientation>::
+computeEquilibrium(int iPop, T rho, const T u[DESCRIPTOR::d], T uSqr) const
+{
+  return boundaryDynamics.computeEquilibrium(iPop, rho, u, uSqr);
+}
+
+template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation>
+void ZouHeDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>::collide (
+  Cell<T,DESCRIPTOR>& cell,
+  LatticeStatistics<T>& statistics )
+{
+  typedef lbHelpers<T,DESCRIPTOR> lbH;
+  typedef DESCRIPTOR L;
+
+  // Along all the commented parts of this code there will be an example based
+  // on the situation where the wall's normal vector if (0,1) and the
+  // numerotation of the velocites are done according to the D2Q9
+  // lattice of the OpenLB library.
+
+  // Find all the missing populations
+  // (directions 3,4,5)
+  std::vector<int> missingIndexes = util::subIndexOutgoing<L,direction,orientation>();
+
+  // Will contain the missing poputations that are not normal to the wall.
+  // (directions 3,5)
+  std::vector<int> missingDiagonalIndexes = missingIndexes;
+  for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
+    int numOfNonNullComp = 0;
+    for (int iDim = 0; iDim < L::d; ++iDim) {
+      numOfNonNullComp += abs(descriptors::c<L>(missingIndexes[iPop],iDim));
+    }
+
+    if (numOfNonNullComp == 1) {
+      missingDiagonalIndexes.erase(missingDiagonalIndexes.begin()+iPop);
+      break;
+    }
+  }
+
+  T rho, u[L::d];
+  T falseRho, falseU[L::d];
+  this->_momenta.computeRhoU(cell, rho, u);
+
+  T uSqr = util::normSqr<T,L::d>(u);
+
+  // The unknown non equilibrium populations are bounced back
+  // (f[3] = feq[3] + fneq[7], f[4] = feq[4] + fneq[8],
+  //  f[5] = feq[5] + fneq[1])
+  for (unsigned iPop = 0; iPop < missingIndexes.size(); ++iPop) {
+    cell[missingIndexes[iPop]] = cell[util::opposite<L>(missingIndexes[iPop])]
+                                 - computeEquilibrium(util::opposite<L>(missingIndexes[iPop]), rho, u, uSqr)
+                                 + computeEquilibrium(missingIndexes[iPop], rho, u, uSqr);
+  }
+
+  // We recompute rho and u in order to have the new momentum and density. Since
+  // the momentum is not conserved from this scheme, we will corect it. By adding
+  // a contribution to the missingDiagonalVelocities.
+  lbH::computeRhoU(cell,falseRho,falseU);
+
+  T diff[L::d];
+  for (int iDim = 0; iDim < L::d; ++iDim) {
+    diff[iDim] = (rho*u[iDim] - falseRho*falseU[iDim])/ (T)missingDiagonalIndexes.size();
+  }
+
+  for (unsigned iPop = 0; iPop < missingDiagonalIndexes.size(); ++iPop) {
+    for (int iDim = 1; iDim < L::d; ++iDim) {
+      cell[missingDiagonalIndexes[iPop]] +=
+        descriptors::c<L>(missingDiagonalIndexes[iPop],(direction+iDim)%L::d) * diff[(direction+iDim)%L::d];
+    }
+  }
+
+  boundaryDynamics.collide(cell, statistics);
+
+  statistics.incrementStats(rho, uSqr);
+}
+
+template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation>
+T ZouHeDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>::getOmega() const
+{
+  return boundaryDynamics.getOmega();
+}
+
+template<typename T, typename DESCRIPTOR, typename Dynamics, int direction, int orientation>
+void ZouHeDynamics<T,DESCRIPTOR,Dynamics,direction,orientation>::setOmega(T omega_)
+{
+  boundaryDynamics.setOmega(omega_);
+}
+
+
+}  // namespace olb
+
+#endif