Initialize at openlb-1-3

1 files changed, 281 insertions, 0 deletions

diff --git a/src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh b/src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh
new file mode 100644
index 0000000..a32d22f
--- /dev/null
+++ b/src/dynamics/navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.hh
@@ -0,0 +1,281 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2008 Orestis Malaspinas, Andrea Parmigiani, 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 NAVIER_STOKES_INTO_ADVECTION_DIFFUSION_MRT_COUPLING_POST_PROCESSOR_3D_HH
+#define NAVIER_STOKES_INTO_ADVECTION_DIFFUSION_MRT_COUPLING_POST_PROCESSOR_3D_HH
+
+#include "latticeDescriptors.h"
+#include "navierStokesAdvectionDiffusionMRTCouplingPostProcessor3D.h"
+#include "core/blockLattice3D.h"
+#include "core/util.h"
+#include "core/finiteDifference3D.h"
+#include "advectionDiffusionForces.hh"
+#include "advectionDiffusionForces.h"
+
+namespace olb {
+
+//=====================================================================================
+//==============  NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D ===========
+//=====================================================================================
+template<typename T, typename DESCRIPTOR>
+NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D<T,DESCRIPTOR>::
+NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
+    T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_,
+    std::vector<SpatiallyExtendedObject3D* > partners_)
+  :  x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
+     gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_),
+     dir(dir_), partners(partners_)
+{
+  // we normalize the direction of force vector
+  T normDir = T();
+  for (unsigned iD = 0; iD < dir.size(); ++iD) {
+    normDir += dir[iD]*dir[iD];
+  }
+  normDir = sqrt(normDir);
+  for (unsigned iD = 0; iD < dir.size(); ++iD) {
+    dir[iD] /= normDir;
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D<T,DESCRIPTOR>::
+processSubDomain(BlockLattice3D<T,DESCRIPTOR>& blockLattice,
+                 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
+{
+  typedef DESCRIPTOR L;
+  enum {x,y,z};
+  enum {
+    velOffset = descriptors::AdvectionDiffusionMRTD3Q7Descriptor::template index<descriptors::VELOCITY>(),
+    forceOffset = DESCRIPTOR::template index<descriptors::FORCE>()
+  };
+
+  BlockLattice3D<T,descriptors::AdvectionDiffusionMRTD3Q7Descriptor> *tPartner =
+    dynamic_cast<BlockLattice3D<T,descriptors::AdvectionDiffusionMRTD3Q7Descriptor> *>(partners[0]);
+
+  int newX0, newX1, newY0, newY1, newZ0, newZ1;
+  if ( util::intersect (
+         x0, x1, y0, y1, z0, z1,
+         x0_, x1_, y0_, y1_, z0_, z1_,
+         newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
+
+    for (int iX=newX0; iX<=newX1; ++iX) {
+      for (int iY=newY0; iY<=newY1; ++iY) {
+        for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
+          //                  Velocity coupling
+          T *u = tPartner->get(iX,iY,iZ).template getFieldPointer<descriptors::VELOCITY>();
+          blockLattice.get(iX,iY,iZ).computeU(u);
+
+          //coupling between the temperature and navier stokes.
+
+          T *force = blockLattice.get(iX,iY,iZ).template getFieldPointer<descriptors::FORCE>();
+          T temperature = tPartner->get(iX,iY,iZ).computeRho();
+          T rho = blockLattice.get(iX,iY,iZ).computeRho();
+          for (unsigned iD = 0; iD < L::d; ++iD) {
+            force[iD] = gravity * rho * (temperature - T0) / deltaTemp * dir[iD];
+          }
+        }
+      }
+    }
+  }
+}
+
+template<typename T, typename DESCRIPTOR>
+void NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D<T,DESCRIPTOR>::
+process(BlockLattice3D<T,DESCRIPTOR>& blockLattice)
+{
+  processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
+}
+
+//=====================================================================================
+//==============  StokesDragCouplingPostProcessor3D ===========
+//=====================================================================================
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+AdvectionDiffusionParticleMRTCouplingPostProcessor3D<T,DESCRIPTOR,ADLattice>::
+AdvectionDiffusionParticleMRTCouplingPostProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int iC_,
+    int offset_, std::vector<SpatiallyExtendedObject3D* > partners_,
+    std::vector<std::reference_wrapper<AdvectionDiffusionForce3D<T, DESCRIPTOR, ADLattice> > > forces_)
+  :  x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_), iC(iC_), offset(offset_), forces(forces_)
+{
+  BlockLattice3D<T,ADLattice> *partnerLattice =
+    dynamic_cast<BlockLattice3D<T,ADLattice> *>(partners_[0]);
+  adCell = &partnerLattice->get(x0,y0,z0);
+  vel = partnerLattice->get(x0,y0,z0)[offset];
+  vel_new = partnerLattice->get(x0,y0,z0)[offset];
+  velXp = partnerLattice->get(x0+1,y0,z0)[offset];
+  velXn = partnerLattice->get(x0-1,y0,z0)[offset];
+  velYp = partnerLattice->get(x0,y0+1,z0)[offset];
+  velYn = partnerLattice->get(x0,y0-1,z0)[offset];
+  velZp = partnerLattice->get(x0,y0,z0+1)[offset];
+  velZn = partnerLattice->get(x0,y0,z0-1)[offset];
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+void AdvectionDiffusionParticleMRTCouplingPostProcessor3D<T,DESCRIPTOR,ADLattice>::
+processSubDomain(BlockLattice3D<T,DESCRIPTOR>& blockLattice,
+                 int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
+{
+  int newX0, newX1, newY0, newY1, newZ0, newZ1;
+
+  int off = (par) ? 3 : 0;
+  int off2 = (par) ? 0 : 3;
+
+  if ( util::intersect (
+         x0, x1, y0, y1, z0, z1,
+         x0_, x1_, y0_, y1_, z0_, z1_,
+         newX0, newX1, newY0, newY1, newZ0, newZ1 ) ) {
+
+    for (int iX=newX0; iX<=newX1; ++iX) {
+      for (int iY=newY0; iY<=newY1; ++iY) {
+        for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
+
+          int latticeR[4] = {iC, iX, iY, iZ};
+          T velGrad[3] = {0.,0.,0.};
+          T forceValue[3] = {0.,0.,0.};
+          T velF[3] = {0.,0.,0.};
+
+          nsCell = &(blockLattice.get(iX,iY,iZ));
+          //.computeU(velF);
+
+          if (forces.begin() != forces.end()) {
+            // calculating upwind Gradient
+            if (vel[0+off]<0.) {
+              velGrad[0] = vel[0+off]*(velXp[0+off]-vel[0+off]);
+              velGrad[1] = vel[0+off]*(velXp[1+off]-vel[1+off]);
+              velGrad[2] = vel[0+off]*(velXp[2+off]-vel[2+off]);
+            } else {
+              velGrad[0] = vel[0+off]*(vel[0+off]-velXn[0+off]);
+              velGrad[1] = vel[0+off]*(vel[1+off]-velXn[1+off]);
+              velGrad[2] = vel[0+off]*(vel[2+off]-velXn[2+off]);
+            }
+            if (vel[1+off]<0.) {
+              velGrad[0] += vel[1+off]*(velYp[0+off]-vel[0+off]);
+              velGrad[1] += vel[1+off]*(velYp[1+off]-vel[1+off]);
+              velGrad[2] += vel[1+off]*(velYp[2+off]-vel[2+off]);
+            } else {
+              velGrad[0] += vel[1+off]*(vel[0+off]-velYn[0+off]);
+              velGrad[1] += vel[1+off]*(vel[1+off]-velYn[1+off]);
+              velGrad[2] += vel[1+off]*(vel[2+off]-velYn[2+off]);
+            }
+            if (vel[2+off]<0.) {
+              velGrad[0] += vel[2+off]*(velZp[0+off]-vel[0+off]);
+              velGrad[1] += vel[2+off]*(velZp[1+off]-vel[1+off]);
+              velGrad[2] += vel[2+off]*(velZp[2+off]-vel[2+off]);
+            } else {
+              velGrad[0] += vel[2+off]*(vel[0+off]-velZn[0+off]);
+              velGrad[1] += vel[2+off]*(vel[1+off]-velZn[1+off]);
+              velGrad[2] += vel[2+off]*(vel[2+off]-velZn[2+off]);
+            }
+
+            for (AdvectionDiffusionForce3D<T, DESCRIPTOR,ADLattice>& f : forces) {
+              f.applyForce(forceValue, nsCell, adCell, &vel[off], latticeR);
+            }
+
+            // compute new particle velocity
+            for (int i=0; i < DESCRIPTOR::d; i++) {
+              vel_new[i+off2] = vel[i+off] + forceValue[i] - velGrad[i];
+            }
+          } else {
+            nsCell->computeU(velF);
+            T scaleFactor = T(1);
+            for (int i = 0; i < DESCRIPTOR::d; i++) {
+              vel_new[i + off2] = velF[i] * scaleFactor;
+            }
+          }
+        }
+      }
+    }
+  }
+  par = !par;
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+void AdvectionDiffusionParticleMRTCouplingPostProcessor3D<T,DESCRIPTOR,ADLattice>::
+process(BlockLattice3D<T,DESCRIPTOR>& blockLattice)
+{
+  processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
+}
+
+// LatticeCouplingGenerator for advectionDiffusion coupling
+
+template<typename T, typename DESCRIPTOR>
+NavierStokesAdvectionDiffusionMRTCouplingGenerator3D<T,DESCRIPTOR>::
+NavierStokesAdvectionDiffusionMRTCouplingGenerator3D(int x0_, int x1_, int y0_, int y1_,int z0_, int z1_,
+    T gravity_, T T0_, T deltaTemp_, std::vector<T> dir_)
+  : LatticeCouplingGenerator3D<T,DESCRIPTOR>(x0_, x1_, y0_, y1_, z0_, z1_),
+    gravity(gravity_), T0(T0_), deltaTemp(deltaTemp_), dir(dir_)
+{ }
+
+template<typename T, typename DESCRIPTOR>
+PostProcessor3D<T,DESCRIPTOR>* NavierStokesAdvectionDiffusionMRTCouplingGenerator3D<T,DESCRIPTOR>::generate (
+  std::vector<SpatiallyExtendedObject3D* > partners) const
+{
+  return new NavierStokesAdvectionDiffusionMRTCouplingPostProcessor3D<T,DESCRIPTOR>(
+           this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, gravity, T0, deltaTemp, dir,partners);
+}
+
+template<typename T, typename DESCRIPTOR>
+LatticeCouplingGenerator3D<T,DESCRIPTOR>* NavierStokesAdvectionDiffusionMRTCouplingGenerator3D<T,DESCRIPTOR>::clone() const
+{
+  return new NavierStokesAdvectionDiffusionMRTCouplingGenerator3D<T,DESCRIPTOR>(*this);
+}
+
+// LatticeCouplingGenerator for one-way advectionDiffusion coupling with Stokes drag
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+AdvectionDiffusionParticleMRTCouplingGenerator3D<T,DESCRIPTOR,ADLattice>::
+AdvectionDiffusionParticleMRTCouplingGenerator3D(int offset_)
+  : LatticeCouplingGenerator3D<T,DESCRIPTOR>(0, 0, 0, 0, 0, 0), offset(offset_)
+{ }
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+PostProcessor3D<T,DESCRIPTOR>* AdvectionDiffusionParticleMRTCouplingGenerator3D<T,DESCRIPTOR,ADLattice>::generate (
+  std::vector<SpatiallyExtendedObject3D* > partners) const
+{
+  return new AdvectionDiffusionParticleMRTCouplingPostProcessor3D<T,DESCRIPTOR,ADLattice>(this->x0,this->x1,this->y0,this->y1,this->z0,this->z1, this->iC, offset, partners, ADforces);
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+LatticeCouplingGenerator3D<T,DESCRIPTOR>* AdvectionDiffusionParticleMRTCouplingGenerator3D<T,DESCRIPTOR,ADLattice>::clone() const
+{
+  return new AdvectionDiffusionParticleMRTCouplingGenerator3D<T,DESCRIPTOR,ADLattice>(*this);
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+void AdvectionDiffusionParticleMRTCouplingGenerator3D<T,DESCRIPTOR,ADLattice>::addForce(
+    AdvectionDiffusionForce3D<T,DESCRIPTOR,ADLattice> &force)
+{
+  ADforces.push_back(force);
+}
+
+
+}  // namespace olb
+
+#endif