diff options
Diffstat (limited to 'src/core/latticeStatistics.hh')
-rw-r--r-- | src/core/latticeStatistics.hh | 344 |
1 files changed, 344 insertions, 0 deletions
diff --git a/src/core/latticeStatistics.hh b/src/core/latticeStatistics.hh new file mode 100644 index 0000000..9fd15b2 --- /dev/null +++ b/src/core/latticeStatistics.hh @@ -0,0 +1,344 @@ +/* This file is part of the OpenLB library + * + * Copyright (C) 2006, 2007, 2008 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. +*/ + +#ifndef LATTICE_STATISTICS_HH +#define LATTICE_STATISTICS_HH + +#include "blockLattice2D.h" +#include "blockLattice3D.h" +#include "latticeStatistics.h" +#include <cmath> +#include <numeric> +#include <limits> +#include "util.h" + +namespace olb { + +////////////////////// Class LatticeStatistics ///////////////// + +template<typename T> +LatticeStatistics<T>::LatticeStatistics() : clout(std::cout,"LatticeStatistics") +{ + initialize(); +} + +template<typename T> +LatticeStatistics<T>::~LatticeStatistics() +{ +} + +template<typename T> +void LatticeStatistics<T>::reset() +{ + // avoid division by zero + if (tmpNumCells == 0) { + for (unsigned iVect=0; iVect<averageVect.size(); ++iVect) { + averageVect[iVect] = T(); + } + for (unsigned iVect=0; iVect<sumVect.size(); ++iVect) { + sumVect[iVect] = T(); + } + for (unsigned iVect=0; iVect<minVect.size(); ++iVect) { + minVect[iVect] = T(); + } + for (unsigned iVect=0; iVect<maxVect.size(); ++iVect) { + maxVect[iVect] = T(); + } + numCells = 0; + firstCall = false; + } else { + // The average density is actually used in the "ConstRhoBgk" model. + // Depending on the simulation setup, it is possible that it has + // a nonsensical value before the simulation is started. For this + // and similar cases, averages are initialized to 1. + if (firstCall) { + for (unsigned iVect=0; iVect<averageVect.size(); ++iVect) { + averageVect[iVect] = (T)1; + } + firstCall = false; + } else { + for (unsigned iVect=0; iVect<averageVect.size(); ++iVect) { + averageVect[iVect] = tmpAv[iVect] / (T)tmpNumCells; + } + } + for (unsigned iVect=0; iVect<sumVect.size(); ++iVect) { + sumVect[iVect] = tmpSum[iVect]; + } + for (unsigned iVect=0; iVect<minVect.size(); ++iVect) { + minVect[iVect] = tmpMin[iVect]; + } + for (unsigned iVect=0; iVect<maxVect.size(); ++iVect) { + maxVect[iVect] = tmpMax[iVect]; + } + averageVect[avEnergy] *= (T)0.5; // energy is 0.5 *uSqr + maxVect[maxU] = sqrt(maxVect[maxU]); // u is sqrt(uSqr) + numCells = tmpNumCells; + } + + for (unsigned iVect=0; iVect<averageVect.size(); ++iVect) { + tmpAv[iVect] = T(); + } + for (unsigned iVect=0; iVect<sumVect.size(); ++iVect) { + tmpSum[iVect] = T(); + } + for (unsigned iVect=0; iVect<minVect.size(); ++iVect) { + tmpMin[iVect] = std::numeric_limits<T>::max(); + } + for (unsigned iVect=0; iVect<maxVect.size(); ++iVect) { + tmpMax[iVect] = std::numeric_limits<T>::min(); + } + + tmpNumCells = 0; +} + +template<typename T> +void LatticeStatistics<T>::reset ( + T average_rho_, T average_energy_, T maxU_, size_t numCells_ ) +{ + averageVect[avRho] = average_rho_; + averageVect[avEnergy] = average_energy_; + maxVect[maxU] = maxU_; + numCells = numCells_; + + tmpAv[avRho] = T(); + tmpAv[avEnergy] = T(); + tmpMax[maxU] = T(); + tmpNumCells = 0; +} + +template<typename T> +void LatticeStatistics<T>::initialize() +{ + tmpAv.resize(2); + averageVect.resize(2); + tmpMax.resize(1); + maxVect.resize(1); + + tmpAv[avRho] = T(); + tmpAv[avEnergy] = T(); + tmpMax[maxU] = T(); + tmpNumCells = 0; + + averageVect[avRho] = (T)1; + averageVect[avEnergy] = T(); + maxVect[maxU] = T(); + + firstCall = true; + + resetTime(); +} + +template<typename T> +int LatticeStatistics<T>::subscribeAverage() +{ + int newSize = tmpAv.size()+1; + tmpAv.resize(newSize); + averageVect.resize(newSize); + return newSize-1; +} + +template<typename T> +int LatticeStatistics<T>::subscribeSum() +{ + int newSize = tmpSum.size()+1; + tmpSum.resize(newSize); + sumVect.resize(newSize); + return newSize-1; +} + +template<typename T> +int LatticeStatistics<T>::subscribeMin() +{ + int newSize = tmpMin.size()+1; + tmpMin.resize(newSize); + minVect.resize(newSize); + return newSize-1; +} + +template<typename T> +int LatticeStatistics<T>::subscribeMax() +{ + int newSize = tmpMax.size()+1; + tmpMax.resize(newSize); + maxVect.resize(newSize); + return newSize-1; +} + +template<typename T> +void LatticeStatistics<T>::incrementStats(T rho, T uSqr) +{ + tmpAv[avRho] += rho; + tmpAv[avEnergy] += uSqr; + if (uSqr > tmpMax[maxU]) { + tmpMax[maxU] = uSqr; + } + ++tmpNumCells; +} + +template<typename T> +void LatticeStatistics<T>::gatherAverage(int whichAverage, T value) +{ + OLB_PRECONDITION( whichAverage < (int) tmpAv.size() ); + tmpAv[whichAverage] += value; +} + +template<typename T> +void LatticeStatistics<T>::gatherSum(int whichSum, T value) +{ + OLB_PRECONDITION( whichSum < (int) tmpSum.size() ); + tmpSum[whichSum] += value; +} + +template<typename T> +void LatticeStatistics<T>::gatherMin(int whichMin, T value) +{ + OLB_PRECONDITION( whichMin < (int) tmpMin.size() ); + if (value < tmpMin[whichMin]) { + tmpMin[whichMin] = value; + } +} + +template<typename T> +void LatticeStatistics<T>::gatherMax(int whichMax, T value) +{ + OLB_PRECONDITION( whichMax < (int) tmpMax.size() ); + if (value > tmpMax[whichMax]) { + tmpMax[whichMax] = value; + } +} + +template<typename T> +void LatticeStatistics<T>::incrementStats() +{ + ++tmpNumCells; +} + +template<typename T> +T LatticeStatistics<T>::getAverageRho() const +{ + return averageVect[avRho]; +} + +template<typename T> +T LatticeStatistics<T>::getAverageEnergy() const +{ + return averageVect[avEnergy]; +} + +template<typename T> +T LatticeStatistics<T>::getMaxU() const +{ + return maxVect[maxU]; +} + +template<typename T> +size_t const& LatticeStatistics<T>::getNumCells() const +{ + return numCells; +} + +template<typename T> +T LatticeStatistics<T>::getAverage(int whichAverage) const +{ + OLB_PRECONDITION( whichAverage < (int) tmpAv.size() ); + return averageVect[whichAverage]; +} + +template<typename T> +T LatticeStatistics<T>::getSum(int whichSum) const +{ + OLB_PRECONDITION( whichSum < (int) tmpSum.size() ); + return sumVect[whichSum]; +} + +template<typename T> +T LatticeStatistics<T>::getMin(int whichMin) const +{ + OLB_PRECONDITION( whichMin < (int) tmpMin.size() ); + return minVect[whichMin]; +} + +template<typename T> +T LatticeStatistics<T>::getMax(int whichMax) const +{ + OLB_PRECONDITION( whichMax < (int) tmpMax.size() ); + return maxVect[whichMax]; +} + +template<typename T> +std::vector<T>& LatticeStatistics<T>::getAverageVect() +{ + return averageVect; +} + +template<typename T> +std::vector<T>& LatticeStatistics<T>::getSumVect() +{ + return sumVect; +} + +template<typename T> +std::vector<T>& LatticeStatistics<T>::getMinVect() +{ + return minVect; +} + +template<typename T> +std::vector<T>& LatticeStatistics<T>::getMaxVect() +{ + return maxVect; +} + +template<typename T> +void LatticeStatistics<T>::incrementTime() +{ + ++latticeTime; +} + +template<typename T> +void LatticeStatistics<T>::resetTime(size_t value) +{ + latticeTime=value; +} + +template<typename T> +size_t LatticeStatistics<T>::getTime() const +{ + return latticeTime; +} + +template<typename T> +void LatticeStatistics<T>::print(int iterationStep, T physicalTime) const +{ + clout + << "step=" << iterationStep << "; " + << "t=" << physicalTime << "; " + << "uMax=" << getMaxU() << "; " + << "avEnergy=" << getAverageEnergy() << "; " + << "avRho=" << getAverageRho() + << std::endl; +} + +} // namespace olb + +#endif |