Initialize at openlb-1-3

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)