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/analytical/analyticalF.hh | 630 +++++++++++++++++++++++++++++++++
1 file changed, 630 insertions(+)
create mode 100644 src/functors/analytical/analyticalF.hh
(limited to 'src/functors/analytical/analyticalF.hh')
diff --git a/src/functors/analytical/analyticalF.hh b/src/functors/analytical/analyticalF.hh
new file mode 100644
index 0000000..83a08a0
--- /dev/null
+++ b/src/functors/analytical/analyticalF.hh
@@ -0,0 +1,630 @@
+/* This file is part of the OpenLB library
+ *
+ * Copyright (C) 2012 Lukas Baron, Tim Dornieden, Mathias J. Krause, Albert Mink
+ * 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 ANALYTICAL_F_HH
+#define ANALYTICAL_F_HH
+
+#include
+#include // for lpnorm
+#include // for random
+#include
+#include
+
+#include "analyticalF.h"
+#include "core/singleton.h"
+#include "communication/mpiManager.h"
+#include "utilities/vectorHelpers.h"
+#include "core/radiativeUnitConverter.h"
+
+
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+namespace olb {
+
+
+//////////////////////////////////1D////////////////////////////////////////////
+template
+AnalyticalConst1D::AnalyticalConst1D(const std::vector& value)
+ : AnalyticalF1D(value.size()), _c(value)
+{
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst1D::AnalyticalConst1D(T value) : AnalyticalF1D(1)
+{
+ _c.push_back(value);
+ this->getName() = "const";
+}
+
+template
+bool AnalyticalConst1D::operator()(T output[], const S x[])
+{
+ for ( unsigned i = 0; i < _c.size(); ++i) {
+ output[i] = _c[i];
+ }
+ return true;
+}
+
+
+template
+AnalyticalLinear1D::AnalyticalLinear1D(T a, T b)
+ : AnalyticalF1D(1), _a(a), _b(b)
+{
+ this->getName() = "linear";
+}
+
+template
+AnalyticalLinear1D::AnalyticalLinear1D(S x0, T v0, S x1, T v1)
+ : AnalyticalF1D(1)
+{
+ if ( util::nearZero(x1-x0) ) {
+ std::cout << "Error: x1-x2=0" << std::endl;
+ } else {
+ _a = ( v1-v0 ) / ( x1-x0 );
+ _b = v0 - _a*x0;
+ }
+ this->getName() = "linear";
+}
+
+template
+bool AnalyticalLinear1D::operator()(T output[], const S x[])
+{
+ output[0]=_a*x[0] + _b;
+ return true;
+}
+
+
+template
+AnalyticalRandom1D::AnalyticalRandom1D() : AnalyticalF1D(1)
+{
+#ifdef PARALLEL_MODE_MPI
+ int nameLen, numProcs, myID;
+ char processorName[MPI_MAX_PROCESSOR_NAME];
+ MPI_Comm_rank(MPI_COMM_WORLD,&myID);
+ MPI_Comm_size(MPI_COMM_WORLD,&numProcs);
+ MPI_Get_processor_name(processorName,&nameLen);
+ srand(time(0)+myID*numProcs + nameLen);
+#else
+ srand(time(nullptr));
+#endif
+ this->getName() = "random";
+}
+
+template
+bool AnalyticalRandom1D::operator()(T output[], const S x[])
+{
+ output[0]=(rand()%RAND_MAX)/(T)RAND_MAX;
+ return true;
+}
+
+
+template
+AnalyticalSquare1D::AnalyticalSquare1D(S cp, S r, T maxi)
+ : AnalyticalF1D(1), _cp(cp), _r(r), _maxi(maxi)
+{
+ this->getName() = "square";
+}
+
+template
+bool AnalyticalSquare1D::operator()(T output[], const S x[])
+{
+ output[0]=_maxi*(1-(x[0]-_cp) * (x[0]-_cp) / (T)_r / (T)_r);
+ return true;
+}
+
+
+
+/////////////////////////////////someOtherFunctors//////////////////////////////
+template
+PolynomialStartScale::PolynomialStartScale(S numTimeSteps, T maxValue)
+ : AnalyticalF1D(1), _numTimeSteps(numTimeSteps), _maxValue(maxValue)
+{
+ this->getName() = "polyStartScale";
+}
+
+template
+bool PolynomialStartScale::operator()(T output[], const S x[])
+{
+ output[0]=(T) x[0] / (T) _numTimeSteps;
+ output[0]=_maxValue * output[0]*output[0]*output[0] * ( 10.0 + output[0] * (6.0*output[0] - 15.0) );
+ return true;
+}
+
+
+template
+SinusStartScale::SinusStartScale(int numTimeSteps, T maxValue)
+ : AnalyticalF1D(1), _numTimeSteps(numTimeSteps), _maxValue(maxValue),
+ _pi(4.0*atan(1.0))
+{
+ this->getName() = "sinusStartScale";
+}
+
+template
+bool SinusStartScale::operator()(T output[], const S x[])
+{
+ output[0]=(_maxValue * (sin(-_pi / 2.0 + (T)x[0] / (T)_numTimeSteps * _pi) + 1.0)) / 2.0;
+ return true;
+}
+
+template
+AnalyticalDiffFD1D::AnalyticalDiffFD1D(AnalyticalF1D& f, T eps) : AnalyticalF1D(1), _f(f), _eps(eps)
+{
+}
+
+template
+bool AnalyticalDiffFD1D::operator() (T output[], const T input[])
+{
+ _f(output,input);
+ T x = output[0];
+ T input2[1];
+ input2[0] = input[0] + _eps;
+ _f(output,input2);
+ output[0] -= x;
+ output[0] /= _eps;
+ return true;
+}
+
+///////////////////////////////////////2D///////////////////////////////////////
+template
+AnalyticalComposed2D::AnalyticalComposed2D( AnalyticalF2D& f0,
+ AnalyticalF2D& f1)
+ : AnalyticalF2D(2), _f0(f0), _f1(f1)
+{
+ this->getName() = "composed";
+}
+
+template
+bool AnalyticalComposed2D::operator()(T output[], const S x[])
+{
+ T tmp1[_f0.getTargetDim()], tmp2[_f1.getTargetDim()];
+ _f0(tmp1,x);
+ _f1(tmp2,x);
+ output[0]=tmp1[0];
+ output[1]=tmp2[0];
+ return true;
+}
+
+
+template
+AnalyticalConst2D::AnalyticalConst2D(const std::vector& value)
+ : AnalyticalF2D(value.size()), _c(value)
+{
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst2D::AnalyticalConst2D(T value) : AnalyticalF2D(1)
+{
+ _c.push_back(value);
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst2D::AnalyticalConst2D(T value0, T value1) : AnalyticalF2D(2)
+{
+ _c.push_back(value0);
+ _c.push_back(value1);
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst2D::AnalyticalConst2D(T value0, T value1, T value2) : AnalyticalF2D(3)
+{
+ _c.push_back(value0);
+ _c.push_back(value1);
+ _c.push_back(value2);
+ this->getName() = "const";
+}
+
+template
+bool AnalyticalConst2D::operator()(T output[], const S x[])
+{
+ for (unsigned i = 0; i < _c.size(); ++i) {
+ output[i] = _c[i];
+ }
+ return true;
+}
+
+
+template
+AnalyticalLinear2D::AnalyticalLinear2D(T a, T b, T c)
+ : AnalyticalF2D(1), _a(a), _b(b), _c(c)
+{
+ this->getName() = "linear";
+}
+
+template
+AnalyticalLinear2D::AnalyticalLinear2D(S x0, S y0, T v0, S x1, S y1,
+ T v1, S x2, S y2, T v2)
+ : AnalyticalF2D(1)
+{
+ this->getName() = "linear";
+ T n2= (x1-x0)*(y2-y0) - (y1-y0)*(x2-x0);
+ if ( util::nearZero(n2) ) {
+ std::cout << "Error function" << std::endl;
+ } else {
+ T n0 = (y1-y0)*(v2-v0) - (v1-v0)*(y2-y0);
+ T n1 = (v1-v0)*(x2-x0) - (x1-x0)*(v2-v0);
+ _a = -n0 / n2;
+ _b = -n1 / n2;
+ _c = (x0*n0 + y0*n1 + v0*n2) / n2;
+ }
+}
+
+template
+bool AnalyticalLinear2D::operator()(T output[], const S x[])
+{
+ output[0]=_a*x[0] + _b*x[1] + _c;
+ return true;
+}
+
+template
+AnalyticalParticleAdsorptionLinear2D::AnalyticalParticleAdsorptionLinear2D(T center[], T radius, T maxValue) : AnalyticalF2D(2), _radius(radius), _maxValue(maxValue)
+{
+ _center[0] = center[0];
+ _center[1] = center[1];
+ this->getName() = "particleAdsorptionLinear2D";
+}
+
+template
+bool AnalyticalParticleAdsorptionLinear2D::operator()(T output[], const S input[])
+{
+ T dist = sqrt((input[0]-_center[0])*(input[0]-_center[0]) + (input[1]-_center[1])*(input[1]-_center[1]));
+
+ if (dist > _radius) {
+ output[0] = T();
+ output[1] = T();
+ return true;
+ } else {
+ output[0] = _maxValue*(T(1) - dist/_radius)*(_center[0]-input[0])/_radius;
+ output[1] = _maxValue*(T(1) - dist/_radius)*(_center[1]-input[1])/_radius;
+ return true;
+ }
+}
+
+
+
+template
+AnalyticalRandom2D::AnalyticalRandom2D() : AnalyticalF2D(1)
+{
+#ifdef PARALLEL_MODE_MPI
+ int nameLen, numProcs, myID;
+ char processorName[MPI_MAX_PROCESSOR_NAME];
+ MPI_Comm_rank(MPI_COMM_WORLD,&myID);
+ MPI_Comm_size(MPI_COMM_WORLD,&numProcs);
+ MPI_Get_processor_name(processorName,&nameLen);
+ srand(time(0)+myID*numProcs + nameLen);
+#else
+ srand(time(nullptr));
+#endif
+ this->getName() = "random";
+}
+
+template
+bool AnalyticalRandom2D::operator()(T output[], const S x[])
+{
+ output[0]=(rand()%RAND_MAX)/(T)RAND_MAX;
+ return true;
+}
+
+template
+ParticleU2D::ParticleU2D(SmoothIndicatorF2D& indicator, UnitConverter const& converter)
+ :AnalyticalF2D(2), _indicator(indicator), _converter(converter)
+{
+ this->getName() = "ParticleU";
+}
+
+template
+bool ParticleU2D::operator()(T output[], const S input[])
+{
+ //two dimensions: u = U + w x r = (Ux, Uy, 0) + (0,0,w) x (X,Y,0) = (Ux, Uy, 0) + (-w*Y, w*X, 0)
+ output[0] = _converter.getLatticeVelocity( _indicator.getVel()[0] - _indicator.getOmega() * (input[1] - _indicator.getPos()[1]) );
+ output[1] = _converter.getLatticeVelocity( _indicator.getVel()[1] + _indicator.getOmega() * (input[0] - _indicator.getPos()[0]) );
+
+ return true;
+}
+
+
+///////////////////////////////////////3D///////////////////////////////////////
+template
+AnalyticalComposed3D::AnalyticalComposed3D(AnalyticalF3D& f0,
+ AnalyticalF3D& f1, AnalyticalF3D& f2)
+ : AnalyticalF3D(3), _f0(f0), _f1(f1), _f2(f2)
+{
+ this->getName() = "composed";
+}
+
+template
+bool AnalyticalComposed3D::operator()(T output[], const S x[])
+{
+ T outputTmp0[_f0.getTargetDim()];
+ _f0(outputTmp0,x);
+ output[0]=outputTmp0[0];
+ T outputTmp1[_f1.getTargetDim()];
+ _f1(outputTmp1,x);
+ output[1]=outputTmp1[0];
+ T outputTmp2[_f2.getTargetDim()];
+ _f2(outputTmp2,x);
+ output[2]=outputTmp2[0];
+ return true;
+}
+
+template
+AnalyticalConst3D::AnalyticalConst3D(const std::vector& value)
+ : AnalyticalF3D(value.size()), _c(value)
+{
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst3D::AnalyticalConst3D(T value) : AnalyticalF3D(1)
+{
+ _c.push_back(value);
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst3D::AnalyticalConst3D(T value0, T value1) : AnalyticalF3D(2)
+{
+ _c.push_back(value0);
+ _c.push_back(value1);
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst3D::AnalyticalConst3D(T value0, T value1, T value2)
+ : AnalyticalF3D(3)
+{
+ _c.push_back(value0);
+ _c.push_back(value1);
+ _c.push_back(value2);
+ this->getName() = "const";
+}
+
+template
+AnalyticalConst3D::AnalyticalConst3D(const Vector& value)
+ : AnalyticalF3D(3)
+{
+ _c.reserve(3);
+ _c.push_back(value[0]);
+ _c.push_back(value[1]);
+ _c.push_back(value[2]);
+ this->getName() = "const";
+}
+
+template
+bool AnalyticalConst3D::operator()(T output[], const S x[])
+{
+ for (unsigned int i = 0; i < _c.size(); ++i) {
+ output[i] = _c[i];
+ }
+ return true;
+}
+
+
+template
+AnalyticalLinear3D::AnalyticalLinear3D(T a, T b, T c, T d)
+ : AnalyticalF3D(1), _a(a), _b(b), _c(c), _d(d)
+{
+ this->getName() = "linear";
+}
+
+template
+AnalyticalLinear3D::AnalyticalLinear3D(S x0, S y0, S z0, T v0, S x1,
+ S y1, S z1, T v1, S x2, S y2, S z2, T v2, S x3, S y3, S z3, T v3)
+ : AnalyticalF3D(1)
+{
+ this->getName() = "linear";
+ T n = ( (y3-y0)*(x1-x0)-(x3-x0)*(y1-y0) ) * ( (x2-x0)*(z1-z0)-(z2-z0)*(x1-x0) )
+ +( (z3-z0)*(x1-x0)-(x3-x0)*(z1-z0) ) * ( (y2-y0)*(x1-x0)-(x2-x0)*(y1-y0) );
+ if ( util::nearZero(n) ) {
+ std::cout << "Error function" << std::endl;
+ } else {
+ T w = ( (y1-y0)*(x3-x0)-(x1-x0)*(y3-y0) ) * ( (v2-v0)-(x2-x0)*(v1-v0) / (x1-x0) )
+ /( (y2-y0)*(x1-x0)-(x2-x0)*(y1-y0) ) + (v3-v0) - (x3-x0)*(v1-v0) / (x1-x0);
+ T zx = (y1-y0)*( (x2-x0)*(z1-z0)-(z2-z0)*(x1-x0) )
+ -(z1-z0)*( (y2-y0)*(x1-x0)-(x2-x0)*(y1-y0) );
+ T rx = (v1-v0)/(x1-x0) - (y1-y0)*(v2-v0) / ( (y2-y0)*(x1-x0)-(x2-x0)*(y1-y0) )
+ +(y1-y0)*(x2-x0)*(v1-v0) / ( (y2-y0)*(x1-x0)*(x1-x0)-(x2-x0)*(y1-y0)*(x1-x0) );
+ T zy = (x1-x0)*( (x2-x0)*(z1-z0)-(z2-z0)*(x1-x0) );
+ T ry = ( (x1-x0)*(v2-v0)-(x2-x0)*(v1-v0) ) / ( (y2-y0)*(x1-x0)-(x2-x0)*(y1-y0) );
+ T zz = (x1-x0)*( (y2-y0)*(x1-x0)-(x2-x0)*(y1-y0) );
+ T h = w/n;
+ _a = rx + zx*h;
+ _b = ry + zy*h;
+ _c = zz*h;
+ _d = v0 - x0*_a - y0*_b - z0*_c;
+ }
+}
+
+template
+bool AnalyticalLinear3D::operator()(T output[], const S x[])
+{
+ output[0]=_a*x[0] + _b*x[1] + _c*x[2] + _d;
+ return true;
+}
+
+
+template
+AnalyticalRandom3D::AnalyticalRandom3D() : AnalyticalF3D(1)
+{
+#ifdef PARALLEL_MODE_MPI
+ int nameLen, numProcs, myID;
+ char processorName[MPI_MAX_PROCESSOR_NAME];
+ MPI_Comm_rank(MPI_COMM_WORLD,&myID);
+ MPI_Comm_size(MPI_COMM_WORLD,&numProcs);
+ MPI_Get_processor_name(processorName,&nameLen);
+ srand(time(0)+myID*numProcs + nameLen);
+#else
+ srand(time(nullptr));
+#endif
+ this->getName() = "random";
+}
+
+template
+bool AnalyticalRandom3D::operator()(T output[], const S x[])
+{
+ output[0]=(rand()%RAND_MAX)/(T)RAND_MAX;
+ return true;
+}
+
+
+template
+AnalyticalScaled3D::AnalyticalScaled3D(AnalyticalF3D& f, T scale)
+ : AnalyticalF3D(f.getTargetDim()), _f(f), _scale(scale)
+{
+ this->getName() = "scaled";
+}
+
+template
+bool AnalyticalScaled3D::operator()(T output[], const S x[])
+{
+ T outputTmp[_f.getTargetDim()];
+ _f(outputTmp,x);
+ for (int iDim = 0; iDim < this->getTargetDim(); ++iDim) {
+ output[iDim] = _scale*outputTmp[iDim];
+ }
+ return true;
+}
+
+
+// see Mink et al. 2016 in Sec.3.1.
+template
+PLSsolution3D::PLSsolution3D(RadiativeUnitConverter const& converter)
+ : AnalyticalF3D(1),
+ _physSigmaEff(std::sqrt( converter.getPhysAbsorption() / converter.getPhysDiffusion() )),
+ _physDiffusionCoefficient(converter.getPhysDiffusion())
+{
+ this->getName() = "PLSsolution3D";
+}
+
+template
+bool PLSsolution3D::operator()(T output[1], const S x[3])
+{
+ double r = std::sqrt( x[0]*x[0] + x[1]*x[1] + x[2]*x[2] );
+ output[0] = 1. / (4.0*M_PI*_physDiffusionCoefficient*r) *std::exp(-_physSigmaEff * r);
+ return true;
+}
+
+
+template
+LightSourceCylindrical3D::LightSourceCylindrical3D(RadiativeUnitConverter const& converter, Vector center)
+ : AnalyticalF3D(1),
+ _physSigmaEff(sqrt( converter.getPhysAbsorption() / converter.getPhysDiffusion() )),
+ _physDiffusionCoefficient(converter.getPhysDiffusion()),
+ _center(center)
+{
+ this->getName() = "LightSourceCylindrical3D";
+}
+
+template
+bool LightSourceCylindrical3D::operator()(T output[1], const S x[3])
+{
+ double r = sqrt( (x[0]-_center[0])*(x[0]-_center[0]) + (x[1]-_center[1])*(x[1]-_center[1]) );
+ if ( util::nearZero(r) ) {
+ std::cout << "Warning: evaluation of \"LightSourceCylindrical3D\" functor close to singularity." << std::endl;
+ }
+ output[0] = 1. / (4.0*M_PI*_physDiffusionCoefficient*r) *exp(-_physSigmaEff * r);
+ return true;
+}
+
+
+template
+Spotlight::Spotlight(Vector position, Vector direction, T falloff)
+ : AnalyticalF3D(1), _position(position), _orientation(direction[0]/norm(direction), direction[1]/norm(direction), direction[2]/norm(direction)), _falloff(falloff)
+{
+ this->getName() = "Spotlight";
+}
+
+template
+bool Spotlight::operator()(T output[1], const S x[3])
+{
+ Vector w(x[0] -_position[0], x[1] -_position[1], x[2] -_position[2]);
+ normalize(w);
+ T cosPhi = w*_orientation;
+ output[0] = 1. / (norm(w)*norm(w)) * std::pow(std::max(0., cosPhi), _falloff);
+ return true;
+}
+
+
+template
+GaussianHill2D::GaussianHill2D(T sigma, Vector x0, T c0)
+ : AnalyticalF2D(1), _sigma(sigma), _x0(x0[0],x0[1]), _c0(c0)
+{
+ this->getName() = "GaussianHill";
+}
+
+template
+bool GaussianHill2D::operator()(T output[1], const S x[2])
+{
+ output[0] = _c0 * exp(- ((x[0]-_x0[0])*(x[0]-_x0[0]) + (x[1]-_x0[1])*(x[1]-_x0[1])) / (2*_sigma*_sigma) );
+ return true;
+}
+
+
+template
+GaussianHillTimeEvolution2D::GaussianHillTimeEvolution2D(T sigma0, T D, T t, Vector x0, Vector u, T c0)
+ : AnalyticalF2D(1), _sigma02(sigma0*sigma0), _D(D), _t(t), _x0(x0[0],x0[1]), _u(u[0],u[1]), _c0(c0)
+{
+ this->getName() = "GaussianHillTimeEvolution";
+}
+
+template
+bool GaussianHillTimeEvolution2D::operator()(T output[1], const S x[2])
+{
+ output[0] = _sigma02 / (_sigma02+2*_D*_t) * _c0 * exp(- ((x[0]-_x0[0]-_u[0]*_t)*(x[0]-_x0[0]-_u[0]*_t) + (x[1]-_x0[1]-_u[1]*_t)*(x[1]-_x0[1]-_u[1]*_t)) / (2*(_sigma02+2*_D*_t) ));
+ return true;
+}
+
+
+template
+ParticleU3D::ParticleU3D(SmoothIndicatorF3D& indicator, UnitConverter const& converter)
+ :AnalyticalF3D(3), _indicator(indicator), _converter(converter)
+{
+ this->getName() = "ParticleU";
+}
+
+template
+bool ParticleU3D::operator()(T output[], const S input[])
+{
+ //three dimensions: u = U + w x r = (Ux, Uy, Uz) + (wx,wy,wz) x (X,Y,Z) = (Ux, Uy, Uz) + (wy*Z-wz*Y, wz*X-wx*Z, wx*Y-wy*X)
+ output[0] = _converter.getLatticeVelocity( _indicator.getVel()[0] +
+ ( _indicator.getOmega()[1]*(input[2] - _indicator.getPos()[2])
+ - _indicator.getOmega()[2]*(input[1] - _indicator.getPos()[1]) ) );
+ output[1] = _converter.getLatticeVelocity( _indicator.getVel()[1] +
+ ( _indicator.getOmega()[2]*(input[0] - _indicator.getPos()[0])
+ - _indicator.getOmega()[0]*(input[2] - _indicator.getPos()[2]) ) );
+ output[2] = _converter.getLatticeVelocity( _indicator.getVel()[2] +
+ ( _indicator.getOmega()[0]*(input[1] - _indicator.getPos()[1])
+ - _indicator.getOmega()[1]*(input[0] - _indicator.getPos()[0]) ) );
+
+ return true;
+}
+
+} // end namespace olb
+
+#endif
--
cgit v1.2.3