diff options
Diffstat (limited to 'src/utilities')
41 files changed, 4369 insertions, 0 deletions
diff --git a/src/utilities/MakeHeader b/src/utilities/MakeHeader new file mode 100644 index 0000000..2086bce --- /dev/null +++ b/src/utilities/MakeHeader @@ -0,0 +1,34 @@ +# This file is part of the OpenLB library +# +# Copyright (C) 2007 Mathias Krause +# 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. + + +generic := + +precompiled := benchmarkUtil \ + timer \ + vectorHelpers \ + arithmetic \ + functorPtr \ + hyperplane2D \ + hyperplane3D \ + hyperplaneLattice2D \ + hyperplaneLattice3D diff --git a/src/utilities/anisoDiscr.h b/src/utilities/anisoDiscr.h new file mode 100644 index 0000000..a509541 --- /dev/null +++ b/src/utilities/anisoDiscr.h @@ -0,0 +1,298 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Julius Woerner, Albert Mink + * 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. +*/ + + +/** \file + * discretization of anisotrphic phasefunction + * DOI 10.1016/j.ijheatmasstransfer.2011.11.009 + * + */ + + +#include <iostream> +#include <vector> +#include <cmath> +#include <stdlib.h> +#include <fstream> +#include <iomanip> +#include <string> +#include <math.h> +#ifdef FEATURE_OPENBLAS +#include <openblas/lapacke.h> +#endif + +#include "utilities/vectorHelpers.h" + +namespace olb { + +/** Function to compute Henyey Greenstein phase funtion + * \param cosTheta cos(theta) of scattering event with net direction theta + * \param g anisotropy factor + */ +double henyeyGreenstein(double cosTheta, double g) +{ + return (1-g*g) / std::pow(1+g*g-2*g*cosTheta ,1.5); +} + +// check for energy conservation +template< int q, typename DESCRIPTOR> +std::array<double,q> testEnergyConservationColumn( const std::array<std::array<double,q>,q>& phi ) +{ + std::array<double,q> discIntegral; + for( int iVec = 0; iVec < q; iVec++ ) { + discIntegral[iVec] = 0; + for( int iSum = 0; iSum < q; iSum++ ) { + discIntegral[iVec] += descriptors::t<double,DESCRIPTOR>(iSum) * phi[iSum][iVec]; + } + } + return discIntegral; +} + +// check for energy conservation +template<int q, typename DESCRIPTOR> +std::array<double,q> testEnergyConservationRow( const std::array<std::array<double,q>,q>& phi ) +{ + std::array<double,q> discIntegral; + for( int iVec = 0; iVec < q; iVec++ ) { + discIntegral[iVec] = 0; + for( int iSum = 0; iSum < q; iSum++ ) { + discIntegral[iVec] += descriptors::t<double,DESCRIPTOR>(iSum) * phi[iVec][iSum]; + } + } + return discIntegral; +} + +// check for anisotropy conservation +template< int q> +std::array<double,q> testAnisotropyConservationColumn( const std::array<std::array<double,q>,q>& phi, + const double weights[q], std::array<std::array<double,q>,q>& cosTheta) +{ + std::array<double,q> discIntegral; + for( int iVec = 0; iVec < q; iVec++ ) { + discIntegral[iVec] = 0; + for( int iSum = 0; iSum < q; iSum++ ) { + discIntegral[iVec] += weights[iSum] * cosTheta[iVec][iSum]* phi[iVec][iSum]; + } + } + return discIntegral; +} + + +/** Computes vector [a, a+stepsize, a+2*stepsize, ..., b] + * \param stepsize + * \param a + * \param b + */ +std::vector<double> linespace( double const stepsize, double const a, double const b ) +{ + std::vector<double> linspace{}; // initalize to empty + if( util::nearZero( a-b ) ) { + return linspace; + } + int const h = (int) ( (b-a) / stepsize); + for (int n = 0; n <= h; n++) { + linspace.push_back( a +double(n)/h * b ); + } + return linspace; +} + +/** + * \tparam DESCRIPTOR a lattice stencil + * \tparam q number of lattice stencils MINUS one, 0th direction not needed + * \param stepsize choose precision + * \param anisotropyFactor also known as g + * \param solution for debugging + * \param phi is anisotropy matrix(sym), element [i][j] probability of scattering from i into j (direction) + * \param breakAfter to limit endless iteration for testing + */ +template<typename DESCRIPTOR> +void computeAnisotropyMatrix( double const stepsize, double const anisotropyFactor, + double solution[(DESCRIPTOR::q-1)*((DESCRIPTOR::q-1)+1)/2], + std::array<std::array<double,DESCRIPTOR::q-1>, DESCRIPTOR::q-1>& phi, int const breakAfter = -1) +{ + using namespace descriptors; + + OstreamManager clout( std::cout, "AnisotropyMatrix" ); + clout << "Call AnistorpiyMatrix()" << std::endl; +#ifdef FEATURE_OPENBLAS + clout << "Compute anisotropy matrix ..." << std::endl; + typedef DESCRIPTOR L; + int const q = L::q-1; + int const mm = 2*q; + int const nn = (q+1)*q/2; + + // qxq matrix 0th row and 0th column are empty + std::vector<std::vector<double>> angleProb; + angleProb.resize(q); + for ( int n = 0; n < q; n++ ) { + angleProb[n].resize(q); + } + double dotProduct; + double normI; + double normJ; + for (int iPop=0; iPop<q; iPop++) { + for (int jPop=0; jPop<q; jPop++) { + // shift by 1 due to notation in DESCRIPTOR/DESCRIPTOR + // exclude 0th direction + dotProduct = c<L>(iPop+1)*c<L>(jPop+1); + normI = std::sqrt( c<L>(iPop+1)*c<L>(iPop+1) ); + normJ = std::sqrt( c<L>(jPop+1)*c<L>(jPop+1) ); + angleProb[iPop][jPop] = dotProduct / (normI*normJ); + } + } + + for (int i=0; i<q; i++) { + for (int j=0; j<q; j++) { + phi[i][j]=1.0; + } + } + + for( int i = 0; i < nn; i++ ) { + solution[i] = 0; + }; + + int iter; + double U[mm][nn]; + double U_array[nn*mm]; + std::vector<double> anisoIterVector = linespace( stepsize, 0, anisotropyFactor ); + + // additional condition only for unit testing + size_t index = 0; + for ( ; index < anisoIterVector.size() && index != (std::size_t)(breakAfter); index++) + { + // wipe matrices and vectors + for (int m = 0; m < mm; m++) { + for (int n = 0; n < nn; n++) { + U[m][n] = 0; + } + } + + // compute matrix U, iter handels symmetry + iter = 0; + for (int m=0; m<q; m++) { + for (int n=m; n<q; n++) { + U[m][iter] = t<double,L>(m+1)*phi[n][m]; + U[n][iter] = t<double,L>(n+1)*phi[m][n]; + U[m+q][iter] = t<double,L>(m+1)*phi[n][m]*angleProb[n][m]; + U[n+q][iter] = t<double,L>(n+1)*phi[m][n]*angleProb[m][n]; + iter++; + } + } + + double sum1; + double sum2; + // compute vector b + for (int n=0; n<q; n++) { + // get sum + sum1 = 0; + sum2 = 0; + for (int k=0; k<q; k++) { + sum1 += t<double,L>(k+1)*phi[k][n]; + sum2 += t<double,L>(k+1)*phi[k][n]*angleProb[k][n]; + } + solution[n] = 1 - sum1; + solution[q+n] = anisoIterVector[index] - sum2; + } + + // transform 2d array to 1d array, column major + for ( int n = 0; n < nn; ++n) { + for ( int m = 0; m < mm; ++m ) { + U_array[n*mm +m] = U[m][n]; + } + } + + //util::print(b, nn, "b"); + + // solve Ax = QRx = R'Q'x = b + // equiv x = Q*R'\x + LAPACKE_dgels( LAPACK_COL_MAJOR, 'N', mm, + nn, 1, U_array, + mm, solution, nn); + //util::print(b, nn, "b after"); + + iter = 0; + for (int m=0; m<q; m++) { + for (int n=m; n<q; n++) { + phi[n][m] = phi[n][m]*(1+solution[iter]); + phi[m][n] = phi[n][m]; + iter++; + } + } + //util::print( testEnergyConservationColumn<q>(phi,L::t), "colum sum elementwise" ); + //util::print( testEnergyConservationRow<q>(phi,L::t), "row sum elementwise" ); + //util::print( testEnergyConservationSec<q>(phi,L::t,angleProb), "second Moment" ); + } + clout << "Terminate after " << index << " iterations" << std::endl; +#endif +} + + +template<typename DESCRIPTOR> +void computeAnisotropyMatrixKimAndLee( double const anisotropyFactor, + std::array<std::array<double,DESCRIPTOR::q>,DESCRIPTOR::q>& phi ) +{ + OstreamManager clout( std::cout, "AnisotropyMatrix_KimAndLee" ); + clout << "Compute anisotropy matrix ..." << std::endl; + typedef DESCRIPTOR L; + int const q = L::q; + + // qxq matrix 0th row and 0th column are empty + std::array<std::array<double,q>, q> cosTheta; + double dotProduct; + double normI; + double normJ; + for (int iPop=0; iPop<q; iPop++) { + for (int jPop=0; jPop<q; jPop++) { + dotProduct = descriptors::c<L>(iPop) * descriptors::c<L>(jPop); + normI = std::sqrt( util::normSqr<int,3>(descriptors::c<L>(iPop)) ); + normJ = std::sqrt( util::normSqr<int,3>(descriptors::c<L>(jPop)) ); + cosTheta[iPop][jPop] = dotProduct /(normI*normJ); + if( util::normSqr<int,3>(descriptors::c<L>(iPop)) == 0 || util::normSqr<int,3>(descriptors::c<L>(jPop)) == 0){ + cosTheta[iPop][jPop] = 0.0; + } + } + } + + std::array<std::array<double,q>, q> phaseFunction; + for (int m=0; m<q; m++) { + for (int n=0; n<q; n++) { + phaseFunction[m][n] = henyeyGreenstein(cosTheta[m][n], anisotropyFactor); + } + } + + for (int m=0; m<q; m++) { + for (int n=0; n<q; n++) { + dotProduct = 0; + for (int i = 0; i < q; ++i) { + dotProduct += phaseFunction[m][i]*descriptors::t<double,L>(i); + } + phi[m][n] = phaseFunction[m][n] / dotProduct; + } + } + //util::print( testEnergyConservationColumn<q>(phi,L::t), "colum sum elementwise" ); + //util::print( testEnergyConservationRow<q>(phi,L::t), "row sum elementwise" ); + //util::print( testAnisotropyConservationColumn<q>(phi,L::t,cosTheta), "anisotropyConservation Moment" ); +} + + +} // namespace olb diff --git a/src/utilities/arithmetic.cpp b/src/utilities/arithmetic.cpp new file mode 100644 index 0000000..c4ea7e0 --- /dev/null +++ b/src/utilities/arithmetic.cpp @@ -0,0 +1,50 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2017 Adrian Kummerlaender + * 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. +*/ + +#include "arithmetic.h" + +namespace olb { + +namespace util { + +template struct minus<int>; +template struct minus<double>; +template struct minus<bool>; + +template struct plus<int>; +template struct plus<double>; +template struct plus<bool>; + +template struct multiplies<int>; +template struct multiplies<double>; +template struct multiplies<bool>; + +template struct divides<int>; +template struct divides<double>; + +template struct power<int>; +template struct power<double>; + +} // end namespace util + +} // end namespace olb diff --git a/src/utilities/arithmetic.h b/src/utilities/arithmetic.h new file mode 100644 index 0000000..53ad061 --- /dev/null +++ b/src/utilities/arithmetic.h @@ -0,0 +1,102 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2017 Adrian Kummerlaender + * 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 UTILITY_ARITHMETIC_H +#define UTILITY_ARITHMETIC_H + +#include <cmath> +#include <functional> + +namespace olb { + +namespace util { + +/// Wrapper of function object std::minus with special handling for bool +/** + * \tparam T Domain of the substraction operation, _without_ operation is + * performed for boolean inputs. + * + * Note that specialization is not required for boolean union (add) and + * intersection (multiply) functors as their behavior is implicitly defined + * by the corresponding arithmetic operation on the integer representation. + **/ +template <typename T> +struct minus { + /// symbol character for functor naming + static const char symbol = '-'; + + constexpr T operator() (const T& lhs, const T& rhs) const + { + return std::minus<T>()(lhs, rhs); + } +}; + +/// Operator specialization for boolean _without_ operation +/** + * i.e. implements what one reasonably expects to happen when substracting + * indicator functors. + **/ +template <> +constexpr bool minus<bool>::operator() (const bool& lhs, const bool& rhs) const +{ + return lhs && !rhs; +} + +/// Wrapper of function object std::plus +template <typename T> +struct plus : public std::plus<T> { + /// symbol character for functor naming + static const char symbol = '+'; +}; + +/// Wrapper of function object std::multiplies +template <typename T> +struct multiplies : public std::multiplies<T> { + /// symbol character for functor naming + static const char symbol = '*'; +}; + +/// Wrapper of function object std::divides +template <typename T> +struct divides : public std::divides<T> { + /// symbol character for functor naming + static const char symbol = '/'; +}; + +/// Power function object +template <typename T> +struct power { + /// symbol character for functor naming + static const char symbol = '^'; + + constexpr T operator() (const T& base, const T& exponent) const + { + return pow(base, exponent); + } +}; + +} // end namespace util + +} // end namespace olb + +#endif diff --git a/src/utilities/benchmarkUtil.cpp b/src/utilities/benchmarkUtil.cpp new file mode 100644 index 0000000..5c9cc1b --- /dev/null +++ b/src/utilities/benchmarkUtil.cpp @@ -0,0 +1,38 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2007 Jonas Latt + * 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. +*/ + +#include "benchmarkUtil.h" +#include "benchmarkUtil.hh" + +namespace olb { + +namespace util { + +template class CircularBuffer<double>; +template class ValueTracer<double>; +template class BisectStepper<double>; +template class Newton1D<double>; + +} + +} diff --git a/src/utilities/benchmarkUtil.h b/src/utilities/benchmarkUtil.h new file mode 100644 index 0000000..3c93b05 --- /dev/null +++ b/src/utilities/benchmarkUtil.h @@ -0,0 +1,144 @@ +#ifndef BENCHMARK_UTIL_H +#define BENCHMARK_UTIL_H + +#include <deque> +#include "io/ostreamManager.h" +#include "functors/analytical/analyticalF.h" + +namespace olb { + +namespace util { + +/// Check time-convergence of a scalar. +/** This class is useful, for example to check convergence of + * the velocity field for the simulation of a stationary flow. + * Convergence is claimed when the standard deviation of the + * monitored value is smaller than epsilon times the average. + * The statistics are taken over a macroscopic time scale of the + * system. + */ +template<typename T> +class ValueTracer { +public: + /// Ctor. + /** \param u The characteristic velocity of the system, for + * computation of the characteristic time scale. + * \param L The characteristic length of the system, for + * computation of the characteristic time scale. + * \param epsilon Precision of the convergence. + */ + ValueTracer(T u, T L, T epsilon); + /// Ctor. + /** \param deltaT corresponds to iteration steps (averaging) + \* \param epsilon allowed derivation of quantity + */ + ValueTracer(int deltaT, T epsilon); + /// Change values of u and L to update characteristic scales of the system. + void resetScale(T u, T L); + /// reinitializes the values + void resetValues(); + /// Get characteristic time scale. + int getDeltaT() const; + /// Feed the object with a new measured scalar. + void takeValue(T val, bool doPrint=false); + /// Test for convergence, with respect to stdDev. + bool hasConverged() const; + /// Test for convergence, with respect to difference between min and max value; + bool convergenceCheck() const; + /// Test for convergence, with respect to difference between min and max value; + bool hasConvergedMinMax() const; + T computeAverage() const; + T computeStdDev(T average) const; + void setEpsilon(T epsilon); +private: + int _deltaT; + T _epsilon; + int _t; + bool _converged; + std::deque<T> _values; + mutable OstreamManager clout; +}; + +/// Propose successive test values of a scalar (e.g. Re) to check stability of a system. +/** At first, the stability limit is explored by constant + * increments/decrements of the scalar, and then, by successive + * bisection. + */ +template<typename T> +class BisectStepper { +public: + /// The only constructor. + /** \param _iniVal Initial guess for the stability limit. + * \param _step Step size at which the value is initially + * incremented/decremented. + */ + BisectStepper(T _iniVal, T _step=0.); + /// Get new value, and indicate if the previous value yielded a stable system or not. + T getVal(bool stable, bool doPrint=false); + /// Test for convergence. + bool hasConverged(T epsilon) const; +private: + T iniVal, currentVal, lowerVal, upperVal; + T step; + enum {first, up, down, bisect} state; + mutable OstreamManager clout; +}; + +/// 1D Newton simple scheme +template<typename T> +class Newton1D { + +protected: + AnalyticalF1D<T,T>& _f; + AnalyticalDiffFD1D<T> _df; + T _yValue; + T _eps; + int _maxIterations; + +public: + Newton1D(AnalyticalF1D<T,T>& f, T yValue = T(), T eps = 1.e-8, int maxIterations = 100); + + T solve(T startValue, bool print=false); +}; + +/// Trapezoidal rule +template<typename T> +class TrapezRuleInt1D { + +protected: + AnalyticalF1D<T,T>& _f; + +public: + TrapezRuleInt1D(AnalyticalF1D<T,T>& f); + + T integrate(T min, T max, int nSteps); +}; + +/** Simple circular buffer to compute average and other quantities + * over pre-defined temporal windows. + * Works with every T supporting += (T or scalar), and /= (scalar) operations, + * including double and Vector<double, size>. + */ +template<typename T> +class CircularBuffer { +public: + CircularBuffer(int size); + /// insert a new entry ed eventually erases the oldest one + void insert(T entry); + /// average over all the entries + T average(); + /// get reference to the last entry for pos=0, the second-to-last for pos=1, and so on + T& get(int pos); + /// return size of the buffer + int getSize(); + +private: + int _size; + std::vector<T> _data; +}; + +} // namespace util + +} // namespace olb + +#endif diff --git a/src/utilities/benchmarkUtil.hh b/src/utilities/benchmarkUtil.hh new file mode 100644 index 0000000..0da96b1 --- /dev/null +++ b/src/utilities/benchmarkUtil.hh @@ -0,0 +1,330 @@ +#ifndef BENCHMARK_UTIL_HH +#define BENCHMARK_UTIL_HH + +#include "benchmarkUtil.h" +#include <iostream> +#include <string> +#include <sstream> +#include <numeric> +#include <algorithm> +#include <cmath> +#include <cassert> + + +namespace olb { + +namespace util { + + +/////////// Class ValueTracer //////////////////////// + +template<typename T> +ValueTracer<T>::ValueTracer(T u, T L, T epsilon) + : _deltaT((int)(L/u/2.)), + _epsilon(epsilon), + _t(0), + _converged(false), + clout(std::cout,"ValueTracer") +{ } + +template<typename T> +ValueTracer<T>::ValueTracer(int deltaT, T epsilon) + : _deltaT(deltaT), + _epsilon(epsilon), + _t(0), + _converged(false), + clout(std::cout,"ValueTracer") +{ } + +template<typename T> +int ValueTracer<T>::getDeltaT() const +{ + return _deltaT; +} + +template<typename T> +void ValueTracer<T>::takeValue(T val, bool doPrint) +{ + _values.push_back(val); + if ((int)_values.size() > abs(_deltaT)) { + _values.erase(_values.begin()); + if (doPrint && _t%_deltaT==0) { + T average = computeAverage(); + T stdDev = computeStdDev(average); + clout << "average=" << average << "; stdDev/average=" << stdDev/average << std::endl; + } + } + ++_t; +} + +template<typename T> +void ValueTracer<T>::resetScale(T u, T L) +{ + _t = _t%_deltaT; + _deltaT = (int) (L/u/2.); + if ( (int)_values.size() > abs(_deltaT) ) { + _values.erase(_values.begin(), _values.begin() + (_values.size()-_deltaT) ); + } +} + +template<typename T> +void ValueTracer<T>::resetValues() +{ + _t = 0; + if ((int)_values.size() > 0) { + _values.erase(_values.begin(), _values.begin() + _values.size() ); + } +} + +template<typename T> +bool ValueTracer<T>::hasConverged() const +{ + if ((int)_values.size() < abs(_deltaT)) { + return false; + } else { + T average = computeAverage(); + T stdDev = computeStdDev(average); + if (!std::isnan(stdDev/average)) { + return fabs(stdDev/average) < _epsilon; + } else { + clout << "simulation diverged." << std::endl; + return true; + } + } +} +template< |