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/utilities/MakeHeader | 34 ++++ src/utilities/anisoDiscr.h | 298 +++++++++++++++++++++++++++++ src/utilities/arithmetic.cpp | 50 +++++ src/utilities/arithmetic.h | 102 ++++++++++ src/utilities/benchmarkUtil.cpp | 38 ++++ src/utilities/benchmarkUtil.h | 144 ++++++++++++++ src/utilities/benchmarkUtil.hh | 330 +++++++++++++++++++++++++++++++++ src/utilities/blockDataReductionMode.h | 45 +++++ src/utilities/blockDataSyncMode.h | 55 ++++++ src/utilities/fraction.h | 73 ++++++++ src/utilities/functorDsl2D.h | 79 ++++++++ src/utilities/functorDsl2D.hh | 112 +++++++++++ src/utilities/functorDsl3D.h | 113 +++++++++++ src/utilities/functorDsl3D.hh | 113 +++++++++++ src/utilities/functorPtr.cpp | 64 +++++++ src/utilities/functorPtr.h | 143 ++++++++++++++ src/utilities/functorPtr.hh | 146 +++++++++++++++ src/utilities/hyperplane2D.cpp | 31 ++++ src/utilities/hyperplane2D.h | 65 +++++++ src/utilities/hyperplane2D.hh | 109 +++++++++++ src/utilities/hyperplane3D.cpp | 31 ++++ src/utilities/hyperplane3D.h | 99 ++++++++++ src/utilities/hyperplane3D.hh | 182 ++++++++++++++++++ src/utilities/hyperplaneLattice2D.cpp | 31 ++++ src/utilities/hyperplaneLattice2D.h | 105 +++++++++++ src/utilities/hyperplaneLattice2D.hh | 199 ++++++++++++++++++++ src/utilities/hyperplaneLattice3D.cpp | 31 ++++ src/utilities/hyperplaneLattice3D.h | 123 ++++++++++++ src/utilities/hyperplaneLattice3D.hh | 300 ++++++++++++++++++++++++++++++ src/utilities/meta.h | 95 ++++++++++ src/utilities/module.mk | 27 +++ src/utilities/namedType.h | 52 ++++++ src/utilities/timer.cpp | 40 ++++ src/utilities/timer.h | 177 ++++++++++++++++++ src/utilities/timer.hh | 266 ++++++++++++++++++++++++++ src/utilities/utilities2D.h | 34 ++++ src/utilities/utilities2D.hh | 33 ++++ src/utilities/utilities3D.h | 35 ++++ src/utilities/utilities3D.hh | 33 ++++ src/utilities/vectorHelpers.cpp | 56 ++++++ src/utilities/vectorHelpers.h | 276 +++++++++++++++++++++++++++ 41 files changed, 4369 insertions(+) create mode 100644 src/utilities/MakeHeader create mode 100644 src/utilities/anisoDiscr.h create mode 100644 src/utilities/arithmetic.cpp create mode 100644 src/utilities/arithmetic.h create mode 100644 src/utilities/benchmarkUtil.cpp create mode 100644 src/utilities/benchmarkUtil.h create mode 100644 src/utilities/benchmarkUtil.hh create mode 100644 src/utilities/blockDataReductionMode.h create mode 100644 src/utilities/blockDataSyncMode.h create mode 100644 src/utilities/fraction.h create mode 100644 src/utilities/functorDsl2D.h create mode 100644 src/utilities/functorDsl2D.hh create mode 100644 src/utilities/functorDsl3D.h create mode 100644 src/utilities/functorDsl3D.hh create mode 100644 src/utilities/functorPtr.cpp create mode 100644 src/utilities/functorPtr.h create mode 100644 src/utilities/functorPtr.hh create mode 100644 src/utilities/hyperplane2D.cpp create mode 100644 src/utilities/hyperplane2D.h create mode 100644 src/utilities/hyperplane2D.hh create mode 100644 src/utilities/hyperplane3D.cpp create mode 100644 src/utilities/hyperplane3D.h create mode 100644 src/utilities/hyperplane3D.hh create mode 100644 src/utilities/hyperplaneLattice2D.cpp create mode 100644 src/utilities/hyperplaneLattice2D.h create mode 100644 src/utilities/hyperplaneLattice2D.hh create mode 100644 src/utilities/hyperplaneLattice3D.cpp create mode 100644 src/utilities/hyperplaneLattice3D.h create mode 100644 src/utilities/hyperplaneLattice3D.hh create mode 100644 src/utilities/meta.h create mode 100644 src/utilities/module.mk create mode 100644 src/utilities/namedType.h create mode 100644 src/utilities/timer.cpp create mode 100644 src/utilities/timer.h create mode 100644 src/utilities/timer.hh create mode 100644 src/utilities/utilities2D.h create mode 100644 src/utilities/utilities2D.hh create mode 100644 src/utilities/utilities3D.h create mode 100644 src/utilities/utilities3D.hh create mode 100644 src/utilities/vectorHelpers.cpp create mode 100644 src/utilities/vectorHelpers.h (limited to 'src/utilities') 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 +# +# +# 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 + * + * + * 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 +#include +#include +#include +#include +#include +#include +#include +#ifdef FEATURE_OPENBLAS +#include +#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 testEnergyConservationColumn( const std::array,q>& phi ) +{ + std::array discIntegral; + for( int iVec = 0; iVec < q; iVec++ ) { + discIntegral[iVec] = 0; + for( int iSum = 0; iSum < q; iSum++ ) { + discIntegral[iVec] += descriptors::t(iSum) * phi[iSum][iVec]; + } + } + return discIntegral; +} + +// check for energy conservation +template +std::array testEnergyConservationRow( const std::array,q>& phi ) +{ + std::array discIntegral; + for( int iVec = 0; iVec < q; iVec++ ) { + discIntegral[iVec] = 0; + for( int iSum = 0; iSum < q; iSum++ ) { + discIntegral[iVec] += descriptors::t(iSum) * phi[iVec][iSum]; + } + } + return discIntegral; +} + +// check for anisotropy conservation +template< int q> +std::array testAnisotropyConservationColumn( const std::array,q>& phi, + const double weights[q], std::array,q>& cosTheta) +{ + std::array 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 linespace( double const stepsize, double const a, double const b ) +{ + std::vector 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 +void computeAnisotropyMatrix( double const stepsize, double const anisotropyFactor, + double solution[(DESCRIPTOR::q-1)*((DESCRIPTOR::q-1)+1)/2], + std::array, 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> 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(iPop+1)*c(jPop+1); + normI = std::sqrt( c(iPop+1)*c(iPop+1) ); + normJ = std::sqrt( c(jPop+1)*c(jPop+1) ); + angleProb[iPop][jPop] = dotProduct / (normI*normJ); + } + } + + for (int i=0; i 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(m+1)*phi[n][m]; + U[n][iter] = t(n+1)*phi[m][n]; + U[m+q][iter] = t(m+1)*phi[n][m]*angleProb[n][m]; + U[n+q][iter] = t(n+1)*phi[m][n]*angleProb[m][n]; + iter++; + } + } + + double sum1; + double sum2; + // compute vector b + for (int n=0; n(k+1)*phi[k][n]; + sum2 += t(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(phi,L::t), "colum sum elementwise" ); + //util::print( testEnergyConservationRow(phi,L::t), "row sum elementwise" ); + //util::print( testEnergyConservationSec(phi,L::t,angleProb), "second Moment" ); + } + clout << "Terminate after " << index << " iterations" << std::endl; +#endif +} + + +template +void computeAnisotropyMatrixKimAndLee( double const anisotropyFactor, + std::array,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, q> cosTheta; + double dotProduct; + double normI; + double normJ; + for (int iPop=0; iPop(iPop) * descriptors::c(jPop); + normI = std::sqrt( util::normSqr(descriptors::c(iPop)) ); + normJ = std::sqrt( util::normSqr(descriptors::c(jPop)) ); + cosTheta[iPop][jPop] = dotProduct /(normI*normJ); + if( util::normSqr(descriptors::c(iPop)) == 0 || util::normSqr(descriptors::c(jPop)) == 0){ + cosTheta[iPop][jPop] = 0.0; + } + } + } + + std::array, q> phaseFunction; + for (int m=0; m(i); + } + phi[m][n] = phaseFunction[m][n] / dotProduct; + } + } + //util::print( testEnergyConservationColumn(phi,L::t), "colum sum elementwise" ); + //util::print( testEnergyConservationRow(phi,L::t), "row sum elementwise" ); + //util::print( testAnisotropyConservationColumn(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 + * + * + * 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; +template struct minus; +template struct minus; + +template struct plus; +template struct plus; +template struct plus; + +template struct multiplies; +template struct multiplies; +template struct multiplies; + +template struct divides; +template struct divides; + +template struct power; +template struct power; + +} // 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 + * + * + * 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 +#include + +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 +struct minus { + /// symbol character for functor naming + static const char symbol = '-'; + + constexpr T operator() (const T& lhs, const T& rhs) const + { + return std::minus()(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::operator() (const bool& lhs, const bool& rhs) const +{ + return lhs && !rhs; +} + +/// Wrapper of function object std::plus +template +struct plus : public std::plus { + /// symbol character for functor naming + static const char symbol = '+'; +}; + +/// Wrapper of function object std::multiplies +template +struct multiplies : public std::multiplies { + /// symbol character for functor naming + static const char symbol = '*'; +}; + +/// Wrapper of function object std::divides +template +struct divides : public std::divides { + /// symbol character for functor naming + static const char symbol = '/'; +}; + +/// Power function object +template +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 + * + * + * 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; +template class ValueTracer; +template class BisectStepper; +template class Newton1D; + +} + +} 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 +#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 +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 _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 +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 +class Newton1D { + +protected: + AnalyticalF1D& _f; + AnalyticalDiffFD1D _df; + T _yValue; + T _eps; + int _maxIterations; + +public: + Newton1D(AnalyticalF1D& f, T yValue = T(), T eps = 1.e-8, int maxIterations = 100); + + T solve(T startValue, bool print=false); +}; + +/// Trapezoidal rule +template +class TrapezRuleInt1D { + +protected: + AnalyticalF1D& _f; + +public: + TrapezRuleInt1D(AnalyticalF1D& 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. + */ +template +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 _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 +#include +#include +#include +#include +#include +#include + + +namespace olb { + +namespace util { + + +/////////// Class ValueTracer //////////////////////// + +template +ValueTracer::ValueTracer(T u, T L, T epsilon) + : _deltaT((int)(L/u/2.)), + _epsilon(epsilon), + _t(0), + _converged(false), + clout(std::cout,"ValueTracer") +{ } + +template +ValueTracer::ValueTracer(int deltaT, T epsilon) + : _deltaT(deltaT), + _epsilon(epsilon), + _t(0), + _converged(false), + clout(std::cout,"ValueTracer") +{ } + +template +int ValueTracer::getDeltaT() const +{ + return _deltaT; +} + +template +void ValueTracer::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 +void ValueTracer::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 +void ValueTracer::resetValues() +{ + _t = 0; + if ((int)_values.size() > 0) { + _values.erase(_values.begin(), _values.begin() + _values.size() ); + } +} + +template +bool ValueTracer::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 +bool ValueTracer::convergenceCheck() 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 false; + } + } +} +template +bool ValueTracer::hasConvergedMinMax() const +{ + if ((int)_values.size() < abs(_deltaT)) { + return false; + } else { + T minEl = *min_element(_values.begin(), _values.end()); + T maxEl = *max_element(_values.begin(), _values.end()); + T average = computeAverage(); + return (maxEl - minEl)/average < _epsilon; + } +} + +template +T ValueTracer::computeAverage() const +{ + return std::accumulate(_values.begin(), _values.end(), 0.) / _values.size(); +} + +template +T ValueTracer::computeStdDev(T average) const +{ + int n = _values.size(); + T sqrDev = 0.; + for (int i=0; i +void ValueTracer::setEpsilon(T epsilon) +{ + _epsilon = epsilon; +} + + + +/////////// Class BisectStepper //////////////////////// + +template +BisectStepper::BisectStepper(T _iniVal, T _step) + : iniVal(_iniVal), step(_step), state(first), clout(std::cout,"BisectStepper") +{ + if (util::nearZero(step)) { + step = iniVal/5.; + } + assert(step>0.); +} + +template +T BisectStepper::getVal(bool stable, bool doPrint) +{ + std::stringstream message; + switch (state) { + case first: + if (stable) { + currentVal = iniVal+step; + state = up; + message << "[" << iniVal << ",infty]"; + } else { + currentVal = iniVal-step; + state = down; + message << "[-infty," << iniVal << "]"; + } + break; + case up: + if (stable) { + message << "[" << currentVal << ",infty]"; + currentVal += step; + } else { + lowerVal = currentVal-step; + upperVal = currentVal; + currentVal = 0.5*(lowerVal+upperVal); + state = bisect; + message << "[" << lowerVal << "," << upperVal << "]"; + } + break; + case down: + if (!stable) { + message << "[-infty," << currentVal << "]"; + currentVal -= step; + } else { + lowerVal = currentVal; + upperVal = currentVal+step; + currentVal = 0.5*(lowerVal+upperVal); + state = bisect; + message << "[" << lowerVal << "," << upperVal << "]"; + } + break; + case bisect: + if (stable) { + lowerVal = currentVal; + } else { + upperVal = currentVal; + } + currentVal = 0.5*(lowerVal+upperVal); + message << "[" << lowerVal << "," << upperVal << "]"; + break; + } + if (doPrint) { + clout << "Value in range " << message.str() << std::endl; + } + return currentVal; +} + +template +bool BisectStepper::hasConverged(T epsilon) const +{ + return (state==bisect) && ((upperVal-lowerVal)/lowerVal < epsilon); +} + +/////////// Class Newton 1D //////////////////////// + +template +Newton1D::Newton1D(AnalyticalF1D& f, T yValue, T eps, int maxIterations) : _f(f), _df(f,eps), _yValue(yValue), _eps(eps), _maxIterations(maxIterations) +{ +} + +template +T Newton1D::solve(T startValue, bool print) +{ + T fValue[1], dfValue[1], x[1]; + x[0] = startValue; + _f(fValue,x); + T eps = fabs(fValue[0] - _yValue); + for (int i=0; i<_maxIterations && eps>_eps; i++) { + _f(fValue,x); + _df(dfValue,x); + eps = fabs(fValue[0] - _yValue); + if (print) { + std::cout << "newtonStep=" << i << "; eps=" << eps << "; x=" << x[0] << "; f(x)="<< fValue[0] << "; df(x)="<< dfValue[0] << std::endl; + } + x[0] -= (fValue[0] - _yValue)/dfValue[0]; + } + return x[0]; +} + +template +TrapezRuleInt1D::TrapezRuleInt1D(AnalyticalF1D& f) : _f(f) +{ +} + +template +T TrapezRuleInt1D::integrate(T min, T max, int nSteps) +{ + T integral[1], tmp[1], min_[1], max_[1], x[1], deltax[0]; + min_[0] = min; + max_[0] = max; + deltax[0] = (max - min) / nSteps; + _f(tmp, min); + integral[0] = 0.5 * tmp[0]; + x[0] = min; + for (int i = 1; i < nSteps; i++) { + x[0] += deltax[0]; + _f(tmp, x); + integral[0] += tmp[0]; + } + + _f(tmp, max); + integral[0] += 0.5*tmp[0]; + integral[0] *= deltax[0]; + + std::cout << "Itnegral=" << integral[0] << std::endl; + + return integral[0]; +} + +/////////// Class CircularBuffer //////////////////////// +template +CircularBuffer::CircularBuffer(int size) + : _size(size) +{} + +template +void CircularBuffer::insert(T entry) +{ + _data.push_back(entry); + if (_data.size() > static_cast(_size) ) { + _data.erase( _data.begin() ); + } +} + +template +T CircularBuffer::average() +{ + T avg = T(); + for (auto i=_data.begin(); i!=_data.end(); i++) { + avg += *i; + } + avg *= 1. / _data.size(); + return avg; +} + +template +T& CircularBuffer::get(int pos) +{ + if (pos < 0) { + return _data.back(); + } + + if (pos >= _size) { + return _data.front(); + } + + return _data[_size - 1 - pos];; +} + +template +int CircularBuffer::getSize() +{ + return _size; +} + + +} // namespace util + +} // namespace olb + +#endif diff --git a/src/utilities/blockDataReductionMode.h b/src/utilities/blockDataReductionMode.h new file mode 100644 index 0000000..38e85b1 --- /dev/null +++ b/src/utilities/blockDataReductionMode.h @@ -0,0 +1,45 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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 BLOCK_DATA_REDUCTION_MODE_H +#define BLOCK_DATA_REDUCTION_MODE_H + +namespace olb { + + +/// Mode of reducing block data from given, possibly higher dimensional data +/** +* Required for optimizing block reduction functors such as BlockReduction3D2D +* if hyperplane is axis-aligned i.e. trivially discretizable. +**/ +enum class BlockDataReductionMode { + /// Interpolate block data at exact physical locations + Analytical, + /// Read block data from discrete lattice locations + Discrete +}; + + +} + +#endif diff --git a/src/utilities/blockDataSyncMode.h b/src/utilities/blockDataSyncMode.h new file mode 100644 index 0000000..9a14a28 --- /dev/null +++ b/src/utilities/blockDataSyncMode.h @@ -0,0 +1,55 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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 BLOCK_DATA_SYNC_MODE_H +#define BLOCK_DATA_SYNC_MODE_H + +namespace olb { + + +/// Mode of synchronizing functor block data between processes +/** +* Required for optimizing functor operations to various usage patterns. +* +* i.e. the convention is for the full domain to be available on every rank +* but this is not ideal in most usage scenarios of reduction functors. +* +* e.g. the primary user of BlockReduction3D2D, BlockGifWriter, only requires +* full data to be available on the rank where its io is performed (rank 0). +* +* SuperLatticeFlux3D only requires rank-local data to be available. Any further +* synchronization would potentially impact performance in this critical area. +**/ +enum class BlockDataSyncMode { + /// default behavior, full block data available on all ranks after update + ReduceAndBcast, + /// optimize for usage in e.g. BlockGifWriter, full data only available on main rank + ReduceOnly, + /// optimize for usage in e.g. SuperLatticeFlux3D, only rank-local data available + None +}; + + +} + +#endif diff --git a/src/utilities/fraction.h b/src/utilities/fraction.h new file mode 100644 index 0000000..2049540 --- /dev/null +++ b/src/utilities/fraction.h @@ -0,0 +1,73 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2019 Adrian Kummerlaender + * 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 UTILITIES_FRACTION_H +#define UTILITIES_FRACTION_H + +#include + +namespace olb { + +namespace utilities { + +/// Floating-point independent fraction type +class Fraction { +private: + const int _numerator; + const int _denominator; + +public: + constexpr Fraction(int num, int denum): + _numerator(num), _denominator(denum) { + if (_denominator == 0) { + throw std::invalid_argument("denominator must not be zero"); + } + } + + constexpr Fraction(int parts[2]): + Fraction(parts[0], parts[1]) { } + + constexpr Fraction(int num): + Fraction(num, 1) { } + + constexpr Fraction(): + Fraction(0) { } + + template + constexpr T as() const + { + return T(_numerator) / T(_denominator); + } + + template + constexpr T inverseAs() const + { + return _numerator != 0 ? T(_denominator) / T(_numerator) : throw std::invalid_argument("inverse of zero is undefined"); + } +}; + +} + +} + +#endif diff --git a/src/utilities/functorDsl2D.h b/src/utilities/functorDsl2D.h new file mode 100644 index 0000000..41ad647 --- /dev/null +++ b/src/utilities/functorDsl2D.h @@ -0,0 +1,79 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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 FUNCTOR_DSL_2D_H +#define FUNCTOR_DSL_2D_H + +#include "functors/lattice/superBaseF2D.h" +#include "functors/analytical/analyticalBaseF.h" + +namespace olb { + +namespace functor_dsl { + +/// Lifts functor reference to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(SuperF2D& f); + +/// Lifts functor pointer to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(SuperF2D* f); + +/// Lifts functor reference to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(AnalyticalF2D& f); + +/// Lifts functor pointer to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(AnalyticalF2D* f); + +/// Returns baseF raised to the power of exponentF +template +std::shared_ptr> pow(std::shared_ptr> baseF, + std::shared_ptr> exponentF); + +/// Returns baseF raised to the power of exponent +template +std::shared_ptr> pow(std::shared_ptr> baseF, + E exponent); + +/// Returns base raised to the power of exponentF +template +std::shared_ptr> pow(B base, + std::shared_ptr> exponentF); + +/// Returns Lp norm for a functor f on the subset described by indicatorF +template +std::shared_ptr> norm(std::shared_ptr> f, + std::shared_ptr> indicatorF); + +/// Returns restriction of a analytical functor f to the lattice sLattice +template +std::shared_ptr> restrict(std::shared_ptr> f, + SuperLattice2D& sLattice); + +} + +} + +#endif diff --git a/src/utilities/functorDsl2D.hh b/src/utilities/functorDsl2D.hh new file mode 100644 index 0000000..7345ef7 --- /dev/null +++ b/src/utilities/functorDsl2D.hh @@ -0,0 +1,112 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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 FUNCTOR_DSL_2D_HH +#define FUNCTOR_DSL_2D_HH + +#include "functorPtr.h" +#include "functors/lattice/superBaseF2D.h" +#include "functors/lattice/indicator/superIndicatorF2D.h" +#include "functors/lattice/integral/superLpNorm2D.h" + +namespace olb { + +namespace functor_dsl { + +template +std::shared_ptr> lift(SuperF2D& f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> lift(SuperF2D* f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> lift(AnalyticalF2D& f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> lift(AnalyticalF2D* f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> pow(std::shared_ptr> baseF, + std::shared_ptr> exponentF) +{ + return std::shared_ptr>( + new SuperCalcPower2D(std::move(baseF), + std::move(exponentF)) + ); +} + +template +std::shared_ptr> pow(std::shared_ptr> baseF, + E exponent) +{ + static_assert(std::is_arithmetic::value, + "Exponent must be an arithmetic value"); + return std::shared_ptr>( + new SuperCalcPower2D(std::move(baseF), exponent)); +} + +template +std::shared_ptr> pow(B base, + std::shared_ptr> exponentF) +{ + static_assert(std::is_arithmetic::value, + "Base must be an arithmetic value"); + return std::shared_ptr>( + new SuperCalcPower2D(base, std::move(exponentF))); +} + +template +std::shared_ptr> norm(std::shared_ptr> f, + std::shared_ptr> indicatorF) +{ + return std::shared_ptr>( + new SuperLpNorm2D(std::move(f), + std::move(indicatorF)) + ); +} + +template +std::shared_ptr> restrict(std::shared_ptr> f, + SuperLattice2D& sLattice) +{ + return std::shared_ptr>( + new SuperLatticeFfromAnalyticalF2D(std::move(f), sLattice)); +} + +} + +} + +#endif diff --git a/src/utilities/functorDsl3D.h b/src/utilities/functorDsl3D.h new file mode 100644 index 0000000..273c90d --- /dev/null +++ b/src/utilities/functorDsl3D.h @@ -0,0 +1,113 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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 FUNCTOR_DSL_3D_H +#define FUNCTOR_DSL_3D_H + +#include "functors/lattice/superBaseF3D.h" +#include "functors/lattice/indicator/superIndicatorF3D.h" +#include "functors/analytical/analyticalBaseF.h" + +namespace olb { + +/// Helper functions for building functors via composition +/** + * Certain types of functors such as error norms (see e.g. SuperRelativeErrorLpNorm3D) + * lend themselves to being expressed as a composition of existing functors. + * + * std::shared_ptr based functor arithmetic offers a memory-safe way to perform + * functor composition. However common operations such as lifting stack-allocated + * functors into the std::shared_ptr context, constructing Lp norms or restricting + * analytical functors on a lattice can become quite verbose. + * + * This is why the `functor_dsl` namespace contains a set of conveniently named + * helper functions that aim to ease writing composed functors. + * + * e.g. a slightly modified version of SuperRelativeErrorLpNorm3D's constructor: + * + * \code{.cpp} + * using namespace functor_dsl; + * + * // decltype(wantedF) == std::shared_ptr> + * // decltype(f) == std::shared_ptr> + * // decltype(indicatorF) == std::shared_ptr> + * + * auto wantedLatticeF = restrict(wantedF, sLattice); + * auto relaticeErrorF = norm

(wantedLatticeF - f, indicatorF) + * / norm

(wantedLatticeF, indicatorF); + * \endcode + * + * Note that any of the functors managed by std::shared_ptr may be returned + * and reused independently. + * + * i.e. it is not necessary to manually assure the lifetime of any functor + * used by the composition. + **/ +namespace functor_dsl { + +/// Lifts functor reference to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(SuperF3D& f); + +/// Lifts functor pointer to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(SuperF3D* f); + +/// Lifts functor reference to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(AnalyticalF3D& f); + +/// Lifts functor pointer to std::shared_ptr functor arithmetic +template +std::shared_ptr> lift(AnalyticalF3D* f); + +/// Returns baseF raised to the power of exponentF +template +std::shared_ptr> pow(std::shared_ptr> baseF, + std::shared_ptr> exponentF); + +/// Returns baseF raised to the power of exponent +template +std::shared_ptr> pow(std::shared_ptr> baseF, + E exponent); + +/// Returns base raised to the power of exponentF +template +std::shared_ptr> pow(B base, + std::shared_ptr> exponentF); + +/// Returns Lp norm for a functor f on the subset described by indicatorF +template +std::shared_ptr> norm(std::shared_ptr> f, + std::shared_ptr> indicatorF); + +/// Returns restriction of a analytical functor f to the lattice sLattice +template +std::shared_ptr> restrict(std::shared_ptr> f, + SuperLattice3D& sLattice); + +} + +} + +#endif diff --git a/src/utilities/functorDsl3D.hh b/src/utilities/functorDsl3D.hh new file mode 100644 index 0000000..85415f7 --- /dev/null +++ b/src/utilities/functorDsl3D.hh @@ -0,0 +1,113 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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 FUNCTOR_DSL_3D_HH +#define FUNCTOR_DSL_3D_HH + +#include "functorPtr.h" +#include "functors/lattice/superBaseF3D.h" +#include "functors/lattice/indicator/superIndicatorF3D.h" +#include "functors/lattice/integral/superLpNorm3D.h" +#include "functors/lattice/reductionF3D.h" + +namespace olb { + +namespace functor_dsl { + +template +std::shared_ptr> lift(SuperF3D& f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> lift(SuperF3D* f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> lift(AnalyticalF3D& f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> lift(AnalyticalF3D* f) +{ + return FunctorPtr>(f).toShared(); +} + +template +std::shared_ptr> pow(std::shared_ptr> baseF, + std::shared_ptr> exponentF) +{ + return std::shared_ptr>( + new SuperCalcPower3D(std::move(baseF), + std::move(exponentF)) + ); +} + +template +std::shared_ptr> pow(std::shared_ptr> baseF, + E exponent) +{ + static_assert(std::is_arithmetic::value, + "Exponent must be an arithmetic value"); + return std::shared_ptr>( + new SuperCalcPower3D(std::move(baseF), exponent)); +} + +template +std::shared_ptr> pow(B base, + std::shared_ptr> exponentF) +{ + static_assert(std::is_arithmetic::value, + "Base must be an arithmetic value"); + return std::shared_ptr>( + new SuperCalcPower3D(base, std::move(exponentF))); +} + +template +std::shared_ptr> norm(std::shared_ptr> f, + std::shared_ptr> indicatorF) +{ + return std::shared_ptr>( + new SuperLpNorm3D(std::move(f), + std::move(indicatorF)) + ); +} + +template +std::shared_ptr> restrict(std::shared_ptr> f, + SuperLattice3D& sLattice) +{ + return std::shared_ptr>( + new SuperLatticeFfromAnalyticalF3D(std::move(f), sLattice)); +} + +} + +} + +#endif diff --git a/src/utilities/functorPtr.cpp b/src/utilities/functorPtr.cpp new file mode 100644 index 0000000..764e2ec --- /dev/null +++ b/src/utilities/functorPtr.cpp @@ -0,0 +1,64 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2018 Adrian Kummerlaender + * 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. +*/ + +#include "functorPtr.h" +#include "functorPtr.hh" + +#include "functors/functors2D.h" +#include "functors/functors3D.h" + +#include "dynamics/latticeDescriptors.h" + +namespace olb { + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>>; + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>>; + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; + +template class FunctorPtr>; +template class FunctorPtr>; +template class FunctorPtr>; + +} diff --git a/src/utilities/functorPtr.h b/src/utilities/functorPtr.h new file mode 100644 index 0000000..438b358 --- /dev/null +++ b/src/utilities/functorPtr.h @@ -0,0 +1,143 @@ +/* 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 + * + * + * 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 FUNCTOR_PTR_H +#define FUNCTOR_PTR_H + +#include +#include + +namespace olb { + +namespace util { + +/// Backport of C++17's std::void_t +template +struct void_t { + using type = void; +}; + +/// Indicates existence of F::identity_functor_type typedef +template +struct has_identity_functor : std::false_type { }; +/// Indicates existence of F::identity_functor_type typedef +template +struct has_identity_functor::type> : std::true_type { }; + +} + +/// Smart pointer for managing the various ways of passing functors around +/** + * There is a rich set of functors that compose other functors by accepting + * them as constructor arguments. e.g. SuperLpNorm3D, SuperPlaneIntegralF3D + * + * Previously all of these functors only accepted other functors by reference + * which prevented e.g. the use case of accepting a std::unique_ptr to a + * indicator freshly created by SuperGeometry3D::getMaterialIndicator. + * + * This class hides the memory management required to support accepting either + * functor references, non-owning pointers or owning pointers in a single argument. + * + * Note that if this class is constructed by reference or raw pointer the caller + * is responsible for assuring that the passed functor is available for the full + * lifetime. This is not the case when constructed by a smart pointer. + **/ +template +class FunctorPtr { +private: + /// Optional functor owner + std::unique_ptr _ownF; + /// Optional shared functor + std::shared_ptr _sharedF; + /// Pointer to the exposed functor + F* _f; + /// True iff FunctorPtr was constructed owning + const bool _owning; + +public: + /// Constructor for transparently accepting a functor reference + FunctorPtr(F& f); + /// Constructor for transparently accepting a non-owning functor pointer + FunctorPtr(F* f); + /// Constructor for transparently accepting a owning functor pointer + /** + * Equivalent to construction using a non-owning functor pointer + **/ + FunctorPtr(const std::unique_ptr& f); + /// Constructor for transparently taking ownership of a owning functor pointer + FunctorPtr(std::unique_ptr&& f); + /// Constructor for transparently sharing ownership of a functor + FunctorPtr(std::shared_ptr f); + /// Constructor for an empty instance + /** + * Equivalent to construction by nullptr. + **/ + FunctorPtr(); + + /// Copy construction is disabled as it is incompatible with unique ownership + FunctorPtr(FunctorPtr&) = delete; + /// Move constructor + FunctorPtr(FunctorPtr&&) = default; + + /// Perfect forwarding functor operator + template + bool operator()(Args... args) + { + return _f->operator()(std::forward(args)...); + } + + /// \returns reference to the exposed functor + typename std::add_lvalue_reference::type operator*() const; + /// Enable pointer-like access to the exposed functor's members + typename std::add_pointer::type operator->() const noexcept; + + /// Indicates whether a functor instance is exposed + /** + * \returns true iff a functor is exposed + * + * This is useful for supporting optional functors in constructor interfaces + * by constructing FunctorPtr using a nullptr + **/ + operator bool() const; + + /// Indicates whether a functor instance is both exposed and owned + /** + * \returns true iff FunctorPtr took ownership of managed functor during construction + * + * It is not enough to check _ownF and _sharedF for non-nullptr contents as _sharedF + * is also used to store F::identity_functor_type instances created by toShared. + **/ + bool isOwning() const; + + /// Lifts managed functor reference for usage in std::shared_ptr arithmetic + template + auto toShared() -> typename std::enable_if< util::has_identity_functor::value, std::shared_ptr>::type; + /// Lifts managed functor reference for usage in std::shared_ptr arithmetic (limited) + template + auto toShared() -> typename std::enable_if::value, std::shared_ptr>::type; + +}; + +} + +#endif diff --git a/src/utilities/functorPtr.hh b/src/utilities/functorPtr.hh new file mode 100644 index 0000000..6d388ef --- /dev/null +++ b/src/utilities/functorPtr.hh @@ -0,0 +1,146 @@ +/* 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 + * + * + * 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 A