summaryrefslogtreecommitdiff
path: root/src/utilities
diff options
context:
space:
mode:
Diffstat (limited to 'src/utilities')
-rw-r--r--src/utilities/MakeHeader34
-rw-r--r--src/utilities/anisoDiscr.h298
-rw-r--r--src/utilities/arithmetic.cpp50
-rw-r--r--src/utilities/arithmetic.h102
-rw-r--r--src/utilities/benchmarkUtil.cpp38
-rw-r--r--src/utilities/benchmarkUtil.h144
-rw-r--r--src/utilities/benchmarkUtil.hh330
-rw-r--r--src/utilities/blockDataReductionMode.h45
-rw-r--r--src/utilities/blockDataSyncMode.h55
-rw-r--r--src/utilities/fraction.h73
-rw-r--r--src/utilities/functorDsl2D.h79
-rw-r--r--src/utilities/functorDsl2D.hh112
-rw-r--r--src/utilities/functorDsl3D.h113
-rw-r--r--src/utilities/functorDsl3D.hh113
-rw-r--r--src/utilities/functorPtr.cpp64
-rw-r--r--src/utilities/functorPtr.h143
-rw-r--r--src/utilities/functorPtr.hh146
-rw-r--r--src/utilities/hyperplane2D.cpp31
-rw-r--r--src/utilities/hyperplane2D.h65
-rw-r--r--src/utilities/hyperplane2D.hh109
-rw-r--r--src/utilities/hyperplane3D.cpp31
-rw-r--r--src/utilities/hyperplane3D.h99
-rw-r--r--src/utilities/hyperplane3D.hh182
-rw-r--r--src/utilities/hyperplaneLattice2D.cpp31
-rw-r--r--src/utilities/hyperplaneLattice2D.h105
-rw-r--r--src/utilities/hyperplaneLattice2D.hh199
-rw-r--r--src/utilities/hyperplaneLattice3D.cpp31
-rw-r--r--src/utilities/hyperplaneLattice3D.h123
-rw-r--r--src/utilities/hyperplaneLattice3D.hh300
-rw-r--r--src/utilities/meta.h95
-rw-r--r--src/utilities/module.mk27
-rw-r--r--src/utilities/namedType.h52
-rw-r--r--src/utilities/timer.cpp40
-rw-r--r--src/utilities/timer.h177
-rw-r--r--src/utilities/timer.hh266
-rw-r--r--src/utilities/utilities2D.h34
-rw-r--r--src/utilities/utilities2D.hh33
-rw-r--r--src/utilities/utilities3D.h35
-rw-r--r--src/utilities/utilities3D.hh33
-rw-r--r--src/utilities/vectorHelpers.cpp56
-rw-r--r--src/utilities/vectorHelpers.h276
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<