src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "correlationFunction.H"
  27
  28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  29
  30 template<class Type>
  31 const char* const
  32     Foam::correlationFunction<Type>::typeName("correlationFunction");
  33
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 template<class Type>
  38 void Foam::correlationFunction<Type>::setTimesAndSizes
  39 (
  40     const label tZeroBufferSize
  41 )
  42 {
  43     sampleSteps_  = ceil(sampleInterval_/mesh_.time().deltaTValue());
  44
  45     sampleInterval_ = sampleSteps_*mesh_.time().deltaTValue();
  46
  47     label bufferLength(ceil(duration_/sampleInterval_));
  48
  49     duration_ = bufferLength*sampleInterval_;
  50
  51     label bufferingInterval(ceil(averagingInterval_/sampleInterval_));
  52
  53     averagingInterval_ = bufferingInterval*sampleInterval_;
  54
  55     label nBuffers(ceil(duration_/averagingInterval_));
  56
  57     this->setSizes
  58     (
  59         nBuffers,
  60         bufferLength,
  61         bufferingInterval
  62     );
  63
  64     tZeroBuffers_ =
  65         Field< Field<Type> >
  66         (
  67             nBuffers,
  68             Field<Type>
  69             (
  70                 tZeroBufferSize,
  71                 pTraits<Type>::zero
  72             )
  73         );
  74 }
  75
  76
  77 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  78
  79 template<class Type>
  80 Foam::correlationFunction<Type>::correlationFunction
  81 (
  82     const polyMesh& mesh,
  83     const dictionary& cfDict,
  84     const label tZeroBufferSize
  85 )
  86 :
  87     bufferedAccumulator<scalar>(),
  88     mesh_(mesh)
  89 {
  90     duration_ = readScalar(cfDict.lookup("duration"));
  91
  92     sampleInterval_ = readScalar(cfDict.lookup("sampleInterval"));
  93
  94     averagingInterval_ = readScalar(cfDict.lookup("averagingInterval"));
  95
  96     setTimesAndSizes(tZeroBufferSize);
  97 }
  98
  99
 100 template<class Type>
 101 Foam::correlationFunction<Type>::correlationFunction
 102 (
 103     const polyMesh& mesh,
 104     const label tZeroBufferSize,
 105     const scalar duration,
 106     const scalar sampleInterval,
 107     const scalar averagingInterval
 108 )
 109 :
 110     bufferedAccumulator<scalar>(),
 111     mesh_(mesh),
 112     duration_(duration),
 113     sampleInterval_(sampleInterval),
 114     averagingInterval_(averagingInterval)
 115 {
 116     setTimesAndSizes(tZeroBufferSize);
 117 }
 118
 119
 120 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 121
 122 template<class Type>
 123 Foam::correlationFunction<Type>::~correlationFunction()
 124 {}
 125
 126
 127 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 128
 129 template<class Type>
 130 void Foam::correlationFunction<Type>::calculateCorrelationFunction
 131 (
 132     const Field<Type>& currentValues
 133 )
 134 {
 135     if (measurandFieldSize() != currentValues.size())
 136     {
 137         FatalErrorIn("correlationFunction<Type>::calculateCorrelationFunction")
 138             << "Trying to supply a Field of length"
 139             << currentValues.size()
 140             << " to calculate the correlation function. "
 141             << "Expecting a Field of length "
 142             << measurandFieldSize() << nl
 143             << abort(FatalError);
 144     }
 145
 146     List<scalar> cFSums(nBuffers(),0.0);
 147
 148     forAll(tZeroBuffers_, tZB)
 149     {
 150         scalar& cFSum = cFSums[tZB];
 151
 152         const Field<Type>& tZeroBuffer = tZeroBuffers_[tZB];
 153
 154         forAll(currentValues, cV)
 155         {
 156             const Type& tZeroBufferValue = tZeroBuffer[cV];
 157
 158             const Type& currentValue = currentValues[cV];
 159
 160             forAll(currentValue, component)
 161             {
 162                 cFSum +=
 163                 (
 164                     tZeroBufferValue[component]*currentValue[component]
 165                 );
 166             }
 167         }
 168
 169         cFSum /= (measurandFieldSize()*currentValues[0].size());
 170     }
 171
 172     label bufferToRefill = addToBuffers(cFSums);
 173
 174     if (bufferToRefill != -1)
 175     {
 176         tZeroBuffers_[bufferToRefill] = currentValues;
 177     }
 178 }
 179
 180
 181 template<class Type>
 182 void Foam::correlationFunction<Type>::calculateCorrelationFunction
 183 (
 184     const Type& currentValue
 185 )
 186 {
 187     if (measurandFieldSize() != 1)
 188     {
 189         FatalErrorIn("correlationFunction<Type>::calculateCorrelationFunction")
 190             << "Trying to supply a single value to calculate the correlation "
 191             << "function.  Expecting a Field of length "
 192             << measurandFieldSize()
 193             << abort(FatalError);
 194     }
 195
 196     calculateCorrelationFunction(Field<Type>(1, currentValue));
 197 }
 198
 199
 200 template<class Type>
 201 Foam::scalar Foam::correlationFunction<Type>::integral() const
 202 {
 203     Field<scalar> averageCF(averaged());
 204
 205     scalar cFIntegral = 0.0;
 206
 207     for (label v = 0; v < averageCF.size() - 1; v++)
 208     {
 209         cFIntegral +=
 210             0.5
 211            *sampleInterval_
 212            *(averageCF[v+1] + averageCF[v]);
 213     }
 214
 215     return cFIntegral;
 216 }
 217
 218
 219 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 220
 221 #   include "correlationFunctionIO.C"
 222
 223 // ************************************************************************* //