Initialize at openlb-1-3

1 files changed, 280 insertions, 0 deletions

diff --git a/src/particles/twoWayCouplings/forwardCouplingModels3D.hh b/src/particles/twoWayCouplings/forwardCouplingModels3D.hh
new file mode 100644
index 0000000..b07d50b
--- /dev/null
+++ b/src/particles/twoWayCouplings/forwardCouplingModels3D.hh
@@ -0,0 +1,280 @@
+/*  Lattice Boltzmann sample, written in C++, using the OpenLB
+ *  library
+ *
+ *  Copyright (C) 2019 Davide Dapelo
+ *  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.
+ */
+
+/* Models for Lagrangian forward-coupling methods -- generic implementation.
+ */
+
+#ifndef LB_FORWARD_COUPLING_MODELS_HH
+#define LB_FORWARD_COUPLING_MODELS_HH
+
+namespace olb {
+
+
+////////////////////// Class LocalBaseCouplingModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+LocalBaseForwardCouplingModel<T,Lattice,Particle>::LocalBaseForwardCouplingModel (
+                     UnitConverter<T, Lattice>& converter,
+                     SuperLattice3D<T, Lattice>& sLattice,
+                     SuperGeometry3D<T>& sGeometry,
+                     DragModel<T, Particle>& dragModel )
+         : _sGeometry(sGeometry),
+           _converter(converter),
+           _sLattice(sLattice),
+           _dragModel(dragModel)
+{
+  this->_interpLatticeDensity = std::make_shared<SuperLatticeInterpDensity3Degree3D<T, Lattice> > (
+                                this->_sLattice, _sGeometry, this->_converter );
+  this->_interpLatticeVelocity = std::make_shared<SuperLatticeInterpPhysVelocity3D<T, Lattice> > (
+                                 this->_sLattice, this->_converter );
+}
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+bool LocalBaseForwardCouplingModel<T,Lattice,Particle>::operator() (Particle<T>* p, int globic)
+{
+  /// Getting the particle and its containing cell's position
+  T physPosP[3] = {T(), T(), T()}; // particle's physical position
+  physPosP[0] = (p->getPos()[0]);
+  physPosP[1] = (p->getPos()[1]);
+  physPosP[2] = (p->getPos()[2]);
+  // particle's dimensionless position, rounded at neighbouring voxel
+  int latticeRoundedPosP[3] = {0, 0, 0};
+  this->_sLattice.getCuboidGeometry().get(globic).getLatticeR (
+                  latticeRoundedPosP, physPosP );
+  // { globic, latticeRoundedP[0, ..., 2] }
+  int globicFull[4] = { globic,
+                        latticeRoundedPosP[0],
+                        latticeRoundedPosP[1],
+                        latticeRoundedPosP[2] };
+
+  // Particle's velocity
+  T physVelP[3] = {T(), T(), T()}; // Physical
+  physVelP[0] = (p->getVel()[0]);
+  physVelP[1] = (p->getVel()[1]);
+  physVelP[2] = (p->getVel()[2]);
+  // Lattice
+  T latticeVelP[3] = {T(), T(), T()}; // particle's dimensionless velocity
+  latticeVelP[0] = this->_converter.getLatticeVelocity(physVelP[0]);
+  latticeVelP[1] = this->_converter.getLatticeVelocity(physVelP[1]);
+  latticeVelP[2] = this->_converter.getLatticeVelocity(physVelP[2]);
+
+  // Lattice's velocity at particle's location
+  T physVelF[3] = {T(), T(), T()}; // Physical
+  this->_interpLatticeVelocity->operator() (physVelF, physPosP, globic);
+  // Lattice
+  T latticeVelF[3] = {T(), T(), T()}; // Lattice's dimensionless velocity at particle's location
+  latticeVelF[0] = this->_converter.getLatticeVelocity(physVelF[0]);
+  latticeVelF[1] = this->_converter.getLatticeVelocity(physVelF[1]);
+  latticeVelF[2] = this->_converter.getLatticeVelocity(physVelF[2]);
+
+  // Computing fluid-particle momentum transfer
+  T gF[3] = {T(), T(), T()}; // force density gF
+  this->_momentumExchange->operator() ( gF, latticeVelF, latticeVelP,
+                           physPosP, latticeRoundedPosP, globic);
+
+  // Computing drag coefficient
+  T Cd = this->_dragModel(p, latticeVelF, latticeVelP, globicFull);
+
+  /// Computing drag force in dimensionless units
+  T latticePRad = p->getRad() / _converter.getConversionFactorLength();
+  T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle
+  latticeForceP[0] = .5 * Cd * M_PI*pow(latticePRad,2) * gF[0] * (latticeVelF[0] - latticeVelP[0]);
+  latticeForceP[1] = .5 * Cd * M_PI*pow(latticePRad,2) * gF[1] * (latticeVelF[1] - latticeVelP[1]);
+  latticeForceP[2] = .5 * Cd * M_PI*pow(latticePRad,2) * gF[2] * (latticeVelF[2] - latticeVelP[2]);
+
+  /// Computing physical drag force
+  std::vector<T> physForceP(3, T()); // physical force acting on the particle
+  physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce();
+  physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce();
+  physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce();
+
+  /// Updating the particle
+  p->setStoreForce(physForceP);
+  return true;
+}
+
+
+////////////////////// Class NaiveForwardCouplingModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+NaiveForwardCouplingModel<T,Lattice,Particle>::NaiveForwardCouplingModel (
+                     UnitConverter<T, Lattice>& converter,
+                     SuperLattice3D<T, Lattice>& sLattice,
+                     SuperGeometry3D<T>& sGeometry,
+                     DragModel<T, Particle>& dragModel )
+         : LocalBaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel)
+{
+  this->_momentumExchange = std::make_shared<NaiveMomentumExchange<T, Lattice> > (
+                            this->_converter, this->_sLattice, this->_interpLatticeDensity );
+}
+
+
+////////////////////// Class LaddForwardCouplingModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+LaddForwardCouplingModel<T,Lattice,Particle>::LaddForwardCouplingModel (
+                     UnitConverter<T, Lattice>& converter,
+                     SuperLattice3D<T, Lattice>& sLattice,
+                     SuperGeometry3D<T>& sGeometry,
+                     DragModel<T, Particle>& dragModel )
+         : LocalBaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel)
+{
+  this->_momentumExchange = std::make_shared<LaddMomentumExchange<T, Lattice> > (
+                            this->_converter, this->_sLattice,
+                            this->_interpLatticeDensity, this->_interpLatticeVelocity );
+}
+
+
+////////////////////// Class NonLocalBaseCouplingModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+NonLocalBaseForwardCouplingModel<T,Lattice,Particle>::NonLocalBaseForwardCouplingModel (
+                     UnitConverter<T, Lattice>& converter,
+                     SuperLattice3D<T, Lattice>& sLattice,
+                     SuperGeometry3D<T>& sGeometry,
+                     DragModel<T, Particle>& dragModel,
+                     SmoothingFunctional<T, Lattice>& smoothingFunctional )
+         : _sGeometry(sGeometry),
+           _converter(converter),
+           _sLattice(sLattice),
+           _dragModel(dragModel),
+           _smoothingFunctional(smoothingFunctional)
+{
+  this->_interpLatticeDensity = std::make_shared<SuperLatticeInterpDensity3Degree3D<T, Lattice> > (
+                                this->_sLattice, _sGeometry, this->_converter );
+  this->_interpLatticeVelocity = std::make_shared<SuperLatticeInterpPhysVelocity3D<T, Lattice> > (
+                                 this->_sLattice, this->_converter );
+}
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+bool NonLocalBaseForwardCouplingModel<T,Lattice,Particle>::operator() (Particle<T>* p, int globic)
+{
+  /// Getting the particle and its containing cell's position
+  T physPosP[3] = {T(), T(), T()}; // particle's physical position
+  physPosP[0] = (p->getPos()[0]);
+  physPosP[1] = (p->getPos()[1]);
+  physPosP[2] = (p->getPos()[2]);
+  //std::cout <<  "globic=" << globic << "   physPosP=(" << physPosP[0] << ", " << physPosP[1] << ", " << physPosP[2] << ")" << std::endl; // <---
+  // particle's dimensionless position, rounded at neighbouring voxel
+  int latticeRoundedPosP[3] = {0, 0, 0};
+  this->_sLattice.getCuboidGeometry().get(globic).getLatticeR (
+           latticeRoundedPosP, physPosP );
+  // { globic, latticeRoundedP[0, ..., 2] }
+  int globicFull[4] = { globic,
+                        latticeRoundedPosP[0],
+                        latticeRoundedPosP[1],
+                        latticeRoundedPosP[2] };
+
+  // Particle's velocity
+  T physVelP[3] = {T(), T(), T()}; // Physical
+  physVelP[0] = (p->getVel()[0]);
+  physVelP[1] = (p->getVel()[1]);
+  physVelP[2] = (p->getVel()[2]);
+  // Lattice
+  T latticeVelP[3] = {T(), T(), T()}; // particle's dimensionless velocity
+  latticeVelP[0] = this->_converter.getLatticeVelocity(physVelP[0]);
+  latticeVelP[1] = this->_converter.getLatticeVelocity(physVelP[1]);
+  latticeVelP[2] = this->_converter.getLatticeVelocity(physVelP[2]);
+    //std::cout <<  "globic=" << globic << "   latticeVelP=(" << latticeVelP[0] << ", " << latticeVelP[1] << ", " << latticeVelP[2] << ")" << std::endl; // <---
+
+  // Update the smoothing functional
+  if ( ! this->_smoothingFunctional.update(physPosP, globic)) {
+    std::cout << "ERROR: no lattice point enclosed in particle's kernel length!" << std::endl;
+    return false;
+  }
+
+  /// Computing dimensionless drag force through smoothed average within the kernel smoothing length
+  T latticeForceP[3] = {T(), T(), T()}; // dimensionless force acting on the particle
+  for (int i=0; i<this->_smoothingFunctional.getSize(); i++) {
+
+    // Position of the iterated voxel
+    int iLatticePosF[3] = {0, 0, 0};
+    this->_smoothingFunctional.getLatticePos(iLatticePosF, i);
+    // Physical
+    T iPhysPosF[3] = {T(), T(), T()};
+    this->_sLattice.getCuboidGeometry().get(globic).getPhysR (
+           iPhysPosF, iLatticePosF );
+
+    // Fluid velocity at the iterated voxel
+    T iPhysVelF[3] = {T(), T(), T()}; // Physical
+    this->_interpLatticeVelocity->operator() (iPhysVelF, iPhysPosF, globic);
+    //std::cout <<  "globic=" << globic << "   iPhysVelF=(" << iPhysVelF[0] << ", " << iPhysVelF[1] << ", " << iPhysVelF[2] << ")" << std::endl; // <---
+    // Lattice
+    T iLatticeVelF[3] = {T(), T(), T()}; // Lattice's dimensionless velocity at particle's location
+    iLatticeVelF[0] = this->_converter.getLatticeVelocity(iPhysVelF[0]);
+    iLatticeVelF[1] = this->_converter.getLatticeVelocity(iPhysVelF[1]);
+    iLatticeVelF[2] = this->_converter.getLatticeVelocity(iPhysVelF[2]);
+
+    // Computing fluid-particle momentum transfer
+    T gF[3] = {T(), T(), T()}; // force density gF
+    this->_momentumExchange->operator() ( gF, iLatticeVelF, latticeVelP,
+                      physPosP, iLatticePosF, globic);
+    //std::cout <<  "globic=" << globic << "   gF=(" << gF[0] << ", " << gF[1] << ", " << gF[2] << ")" << std::endl; // <---
+
+    // Computing drag coefficient
+    T Cd = this->_dragModel(p, iLatticeVelF, latticeVelP, globicFull);
+    //std::cout <<  "globic=" << globic << "   Cd=" << Cd << std::endl; //<---
+
+    /// Computing drag force in dimensionless units
+    T latticePRad = p->getRad() / _converter.getConversionFactorLength();
+    latticeForceP[0] += .5 * Cd * M_PI*pow(latticePRad,2) * gF[0] * (iLatticeVelF[0] - latticeVelP[0])
+                      * this->_smoothingFunctional.getWeight(i);
+    latticeForceP[1] += .5 * Cd * M_PI*pow(latticePRad,2) * gF[1] * (iLatticeVelF[1] - latticeVelP[1])
+                      * this->_smoothingFunctional.getWeight(i);
+    latticeForceP[2] += .5 * Cd * M_PI*pow(latticePRad,2) * gF[2] * (iLatticeVelF[2] - latticeVelP[2])
+                      * this->_smoothingFunctional.getWeight(i);
+  }
+
+  /// Computing physical drag force
+  std::vector<T> physForceP(3, T()); // physical force acting on the particle
+  physForceP[0] = latticeForceP[0] * this->_converter.getConversionFactorForce();
+  physForceP[1] = latticeForceP[1] * this->_converter.getConversionFactorForce();
+  physForceP[2] = latticeForceP[2] * this->_converter.getConversionFactorForce();
+
+  /// Updating the particle
+  p->setStoreForce(physForceP);
+  //std::cout <<  "globic=" << globic << "   physForceP=(" << physForceP[0] << ", " << physForceP[1] << ", " << physForceP[2] << ")" << std::endl; // <---
+  return true;
+}
+
+
+////////////////////// Class NaiveNonLocalForwardCouplingModel ////////////////////////
+
+template<typename T, typename Lattice, template<typename V> class Particle>
+NaiveNonLocalForwardCouplingModel<T,Lattice,Particle>::NaiveNonLocalForwardCouplingModel (
+                     UnitConverter<T, Lattice>& converter,
+                     SuperLattice3D<T, Lattice>& sLattice,
+                     SuperGeometry3D<T>& sGeometry,
+                     DragModel<T, Particle>& dragModel,
+                     SmoothingFunctional<T, Lattice>& smoothingFunctional )
+         : NonLocalBaseForwardCouplingModel<T,Lattice,Particle>(converter, sLattice, sGeometry, dragModel, smoothingFunctional)
+{
+  this->_momentumExchange = std::make_shared<NaiveMomentumExchange<T, Lattice> > (
+                            this->_converter, this->_sLattice, this->_interpLatticeDensity );
+}
+
+
+}
+
+#endif