Initialize at openlb-1-3

1 files changed, 160 insertions, 0 deletions

diff --git a/src/functors/lattice/latticeFrameChangeF3D.hh b/src/functors/lattice/latticeFrameChangeF3D.hh
new file mode 100644
index 0000000..709a76b
--- /dev/null
+++ b/src/functors/lattice/latticeFrameChangeF3D.hh
@@ -0,0 +1,160 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2013, 2015 Gilles Zahnd, Mathias J. Krause
+ *  Marie-Luise Maier
+ *  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 LATTICE_FRAME_CHANGE_F_3D_HH
+#define LATTICE_FRAME_CHANGE_F_3D_HH
+
+#include<cmath>
+
+#include "latticeFrameChangeF3D.h"
+#include "functors/analytical/frameChangeF2D.h"
+#include "core/superLattice3D.h"
+#include "dynamics/lbHelpers.h"  // for computation of lattice rho and velocity
+#include "utilities/vectorHelpers.h"  // for normalize
+#include "geometry/superGeometry3D.h"
+
+namespace olb {
+
+
+template <typename T, typename DESCRIPTOR>
+RotatingForceField3D<T,DESCRIPTOR>::RotatingForceField3D
+(SuperLattice3D<T,DESCRIPTOR>& sLattice_, SuperGeometry3D<T>& superGeometry_,
+ const UnitConverter<T,DESCRIPTOR>& converter_, std::vector<T> axisPoint_,
+ std::vector<T> axisDirection_, T w_, bool centrifugeForceOn_,
+ bool coriolisForceOn_)
+  : SuperLatticeF3D<T,DESCRIPTOR>(sLattice_,3), sg(superGeometry_),
+    converter(converter_), axisPoint(axisPoint_), axisDirection(axisDirection_),
+    w(w_), centrifugeForceOn(centrifugeForceOn_), coriolisForceOn(coriolisForceOn_),
+    velocity(sLattice_,converter_), rho(sLattice_)
+{
+  this->getName() = "rotatingForce";
+}
+
+
+template <typename T, typename DESCRIPTOR>
+bool RotatingForceField3D<T,DESCRIPTOR>::operator()(T output[], const int x[])
+{
+  std::vector<T> F_centri(3,0);
+  std::vector<T> F_coriolis(3,0);
+
+  if ( this->_sLattice.getLoadBalancer().rank(x[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]),&(x[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] = w*w*(physR[0]-axisPoint[0]-scalar*axisDirection[0]);
+      F_centri[1] = w*w*(physR[1]-axisPoint[1]-scalar*axisDirection[1]);
+      F_centri[2] = w*w*(physR[2]-axisPoint[2]-scalar*axisDirection[2]);
+    }
+    if (coriolisForceOn) {
+      T _vel[3];
+      (velocity)(_vel,x);
+      F_coriolis[0] = -2*w*(axisDirection[1]*_vel[2]-axisDirection[2]*_vel[1]);
+      F_coriolis[1] = -2*w*(axisDirection[2]*_vel[0]-axisDirection[0]*_vel[2]);
+      F_coriolis[2] = -2*w*(axisDirection[0]*_vel[1]-axisDirection[1]*_vel[0]);
+    }
+    // return latticeForce
+    output[0] = (F_coriolis[0]+F_centri[0]) * converter.getConversionFactorTime() / converter.getConversionFactorVelocity();
+    output[1] = (F_coriolis[1]+F_centri[1]) * converter.getConversionFactorTime() / converter.getConversionFactorVelocity();
+    output[2] = (F_coriolis[2]+F_centri[2]) * converter.getConversionFactorTime() / converter.getConversionFactorVelocity();
+  }
+  return true;
+}
+
+
+template <typename T, typename DESCRIPTOR>
+HarmonicOscillatingRotatingForceField3D<T,DESCRIPTOR>::HarmonicOscillatingRotatingForceField3D
+(SuperLattice3D<T,DESCRIPTOR>& sLattice_, SuperGeometry3D<T>& superGeometry_,
+ const UnitConverter<T,DESCRIPTOR>& converter_, std::vector<T> axisPoint_,
+ std::vector<T> axisDirection_, T amplitude_, T frequency_)
+  : SuperLatticeF3D<T,DESCRIPTOR>(sLattice_,3), sg(superGeometry_),
+    converter(converter_), axisPoint(axisPoint_), axisDirection(axisDirection_),
+    amplitude(amplitude_), resonanceFrequency(2.*4.*std::atan(1.0)*frequency_), w(0.0), dwdt(0.0),
+    velocity(sLattice_,converter_)
+{
+  this->getName() = "harmonicOscillatingrotatingForce";
+}
+
+template <typename T, typename DESCRIPTOR>
+void HarmonicOscillatingRotatingForceField3D<T,DESCRIPTOR>::updateTimeStep(int iT) {
+  w = resonanceFrequency * amplitude * cos(resonanceFrequency*converter.getPhysTime(iT));
+  dwdt = -resonanceFrequency * resonanceFrequency * amplitude * sin(resonanceFrequency*converter.getPhysTime(iT));
+}
+
+
+template <typename T, typename DESCRIPTOR>
+bool HarmonicOscillatingRotatingForceField3D<T,DESCRIPTOR>::operator()(T output[], const int x[])
+{
+
+
+  std::vector<T> F_centri(3,0);
+  std::vector<T> F_coriolis(3,0);
+  std::vector<T> F_euler(3,0);
+
+  if ( this->_sLattice.getLoadBalancer().rank(x[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]),&(x[0]));
+
+    T scalar =  (physR[0]-axisPoint[0])*axisDirection[0]
+                +(physR[1]-axisPoint[1])*axisDirection[1]
+                +(physR[2]-axisPoint[2])*axisDirection[2];
+
+    T r[3];
+    r[0] = physR[0]-axisPoint[0]-scalar*axisDirection[0];
+    r[1] = physR[1]-axisPoint[1]-scalar*axisDirection[1];
+    r[2] = physR[2]-axisPoint[2]-scalar*axisDirection[2];
+
+    F_centri[0] = w*w*(r[0]);
+    F_centri[1] = w*w*(r[1]);
+    F_centri[2] = w*w*(r[2]);
+
+    T _vel[3];
+    (velocity)(_vel,x);
+    F_coriolis[0] = -2*w*(axisDirection[1]*_vel[2]-axisDirection[2]*_vel[1]);
+    F_coriolis[1] = -2*w*(axisDirection[2]*_vel[0]-axisDirection[0]*_vel[2]);
+    F_coriolis[2] = -2*w*(axisDirection[0]*_vel[1]-axisDirection[1]*_vel[0]);
+
+    F_euler[0] = -dwdt*(axisDirection[1]*r[2]-axisDirection[2]*r[1]);
+    F_euler[1] = -dwdt*(axisDirection[2]*r[0]-axisDirection[0]*r[2]);
+    F_euler[2] = -dwdt*(axisDirection[0]*r[1]-axisDirection[1]*r[0]);
+
+
+    // return latticeForce
+    output[0] = (F_coriolis[0]+F_centri[0]+F_euler[0]) * converter.getConversionFactorTime() / converter.getConversionFactorVelocity();
+    output[1] = (F_coriolis[1]+F_centri[1]+F_euler[1]) * converter.getConversionFactorTime() / converter.getConversionFactorVelocity();
+    output[2] = (F_coriolis[2]+F_centri[2]+F_euler[2]) * converter.getConversionFactorTime() / converter.getConversionFactorVelocity();
+  }
+  return true;
+}
+
+
+} // end namespace olb
+
+#endif