From 94d3e79a8617f88dc0219cfdeedfa3147833719d Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 24 Jun 2019 14:43:36 +0200
Subject: Initialize at openlb-1-3
---
src/functors/lattice/latticeFrameChangeF3D.hh | 160 ++++++++++++++++++++++++++
1 file changed, 160 insertions(+)
create mode 100644 src/functors/lattice/latticeFrameChangeF3D.hh
(limited to 'src/functors/lattice/latticeFrameChangeF3D.hh')
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
+ *
+ *
+ * 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
+
+#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
+RotatingForceField3D::RotatingForceField3D
+(SuperLattice3D& sLattice_, SuperGeometry3D& superGeometry_,
+ const UnitConverter& converter_, std::vector axisPoint_,
+ std::vector axisDirection_, T w_, bool centrifugeForceOn_,
+ bool coriolisForceOn_)
+ : SuperLatticeF3D(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
+bool RotatingForceField3D::operator()(T output[], const int x[])
+{
+ std::vector F_centri(3,0);
+ std::vector 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 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
+HarmonicOscillatingRotatingForceField3D::HarmonicOscillatingRotatingForceField3D
+(SuperLattice3D& sLattice_, SuperGeometry3D& superGeometry_,
+ const UnitConverter& converter_, std::vector axisPoint_,
+ std::vector axisDirection_, T amplitude_, T frequency_)
+ : SuperLatticeF3D(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
+void HarmonicOscillatingRotatingForceField3D::updateTimeStep(int iT) {
+ w = resonanceFrequency * amplitude * cos(resonanceFrequency*converter.getPhysTime(iT));
+ dwdt = -resonanceFrequency * resonanceFrequency * amplitude * sin(resonanceFrequency*converter.getPhysTime(iT));
+}
+
+
+template
+bool HarmonicOscillatingRotatingForceField3D::operator()(T output[], const int x[])
+{
+
+
+ std::vector F_centri(3,0);
+ std::vector F_coriolis(3,0);
+ std::vector 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 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
--
cgit v1.2.3