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