Initialize at openlb-1-3

1 files changed, 130 insertions, 0 deletions

diff --git a/src/dynamics/advectionDiffusionForces.hh b/src/dynamics/advectionDiffusionForces.hh
new file mode 100644
index 0000000..aed1a78
--- /dev/null
+++ b/src/dynamics/advectionDiffusionForces.hh
@@ -0,0 +1,130 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2016 Robin Trunk
+ *  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 ADVECTION_DIFFUSION_FORCES_HH
+#define ADVECTION_DIFFUSION_FORCES_HH
+
+#include "advectionDiffusionForces.h"
+
+namespace olb {
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+AdvDiffDragForce3D<T,DESCRIPTOR,ADLattice>::AdvDiffDragForce3D(UnitConverter<T,DESCRIPTOR> const& converter_, T St_)
+{
+  initArg = 8;
+  dragCoeff = (converter_.getCharPhysVelocity()*converter_.getConversionFactorTime()) / (St_ * converter_.getCharPhysLength());
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+AdvDiffDragForce3D<T,DESCRIPTOR,ADLattice>::AdvDiffDragForce3D(UnitConverter<T,DESCRIPTOR> const& converter_, T pRadius_, T pRho_)
+{
+  initArg = 8;
+  dragCoeff = (9.*converter_.getPhysViscosity()*converter_.getPhysDensity()*converter_.getConversionFactorTime()) / (2.*pRho_*pRadius_*pRadius_);
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+void AdvDiffDragForce3D<T,DESCRIPTOR,ADLattice>::applyForce(T force[], Cell<T,DESCRIPTOR> *nsCell, Cell<T,ADLattice> *adCell, T vel[], int latticeR[])
+{
+  T velF[3] = {0.,0.,0.};
+  nsCell->computeU(velF);
+  for (int i=0; i < DESCRIPTOR::d; i++) {
+    force[i] += dragCoeff*(velF[i]-vel[i]);
+  }
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+AdvDiffRotatingForce3D<T,DESCRIPTOR,ADLattice>::AdvDiffRotatingForce3D(SuperGeometry3D<T>& superGeometry_,
+    const UnitConverter<T,DESCRIPTOR>& converter_, std::vector<T> axisPoint_, std::vector<T> axisDirection_,
+    T w_, T* frac_, bool centrifugeForceOn_, bool coriolisForceOn_) :
+  sg(superGeometry_), axisPoint(axisPoint_), axisDirection(axisDirection_),
+  w(w_), frac(frac_), centrifugeForceOn(centrifugeForceOn_), coriolisForceOn(coriolisForceOn_)
+{
+  invMassLessForce = converter_.getConversionFactorTime() * converter_.getConversionFactorTime() / converter_.getConversionFactorLength();
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+void AdvDiffRotatingForce3D<T,DESCRIPTOR,ADLattice>::applyForce(T force[], Cell<T,DESCRIPTOR> *nsCell, Cell<T,ADLattice> *adCell, T vel[], int latticeR[])
+{
+  std::vector<T> F_centri(3,0);
+  std::vector<T> F_coriolis(3,0);
+  T wf = w*(*frac);
+//  if ( this->_sLattice.getLoadBalancer().rank(latticeR[0]) == singleton::mpi().getRank() ) {
+  // local coords are given, fetch local cell and compute value(s)
+  std::vector<T> physR(3,T());
+  this->sg.getCuboidGeometry().getPhysR(&(physR[0]),&(latticeR[0]));
+
+  T scalar =  (physR[0]-axisPoint[0])*axisDirection[0]
+              +(physR[1]-axisPoint[1])*axisDirection[1]
+              +(physR[2]-axisPoint[2])*axisDirection[2];
+
+  if (centrifugeForceOn) {
+    F_centri[0] = wf*wf*(physR[0]-axisPoint[0]-scalar*axisDirection[0]);
+    F_centri[1] = wf*wf*(physR[1]-axisPoint[1]-scalar*axisDirection[1]);
+    F_centri[2] = wf*wf*(physR[2]-axisPoint[2]-scalar*axisDirection[2]);
+  }
+  if (coriolisForceOn) {
+    F_coriolis[0] = -2*wf*(axisDirection[1]*vel[2]-axisDirection[2]*vel[1]);
+    F_coriolis[1] = -2*wf*(axisDirection[2]*vel[0]-axisDirection[0]*vel[2]);
+    F_coriolis[2] = -2*wf*(axisDirection[0]*vel[1]-axisDirection[1]*vel[0]);
+  }
+  force[0] += (F_coriolis[0]+F_centri[0])*invMassLessForce;
+  force[1] += (F_coriolis[1]+F_centri[1])*invMassLessForce;
+  force[2] += (F_coriolis[2]+F_centri[2])*invMassLessForce;
+//  }
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+AdvDiffMagneticWireForce3D<T,DESCRIPTOR,ADLattice>::AdvDiffMagneticWireForce3D(SuperGeometry3D<T>& superGeometry_, UnitConverter<T,DESCRIPTOR> const& converter_, T pMass, AnalyticalF3D<T, T>& getMagForce) : sg(superGeometry_), _getMagForce(getMagForce)
+{
+  initArg = 8;
+  _pMass = converter_.getConversionFactorTime() / pMass;
+  _conversionVelocity = converter_.getConversionFactorVelocity();
+}
+
+template<typename T, typename DESCRIPTOR,
+typename ADLattice>
+void AdvDiffMagneticWireForce3D<T,DESCRIPTOR,ADLattice>::applyForce(T force[], Cell<T,DESCRIPTOR> *nsCell, Cell<T,ADLattice> *adCell, T vel[], int latticeR[])
+{
+  std::vector<T> physR(3,T());
+  this->sg.getCuboidGeometry().getPhysR(&(physR[0]),&(latticeR[0]));
+  T pos[3] = { T(), T(), T() };
+  pos[0] = physR[0];
+  pos[1] = physR[1];
+  pos[2] = physR[2];
+  T forceHelp[3] = { T(), T(), T() };
+  _getMagForce(forceHelp, pos);
+  for (int i=0; i < DESCRIPTOR::d; i++) {
+    force[i] += forceHelp[i] * _pMass / _conversionVelocity;
+//    std::cout << "----->>>>> force " << forceHelp[i] << std::endl;
+  }
+}
+
+
+}
+#endif