src/randomProcesses/fft/kShellIntegration.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "kShellIntegration.H"
  28 #include "mathematicalConstants.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34
  35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  36
  37 graph kShellIntegration
  38 (
  39     const complexVectorField& Ek,
  40     const Kmesh& K
  41 )
  42 {
  43     // evaluate the radial component of the spectra as an average
  44     // over the shells of thickness dk
  45
  46     graph kShellMeanEk = kShellMean(Ek, K);
  47     const scalarField& x = kShellMeanEk.x();
  48     scalarField& y = *kShellMeanEk.begin()();
  49
  50     // now multiply by 4pi k^2 (the volume of each shell) to get the
  51     // spectra E(k). int E(k) dk is now the total energy in a box
  52     // of side 2pi
  53
  54     y *= sqr(x)*4.0*mathematicalConstant::pi;
  55
  56     // now scale this to get the energy in a box of side l0
  57
  58     scalar l0(K.sizeOfBox()[0]*(scalar(K.nn()[0])/(scalar(K.nn()[0])-1.0)));
  59     scalar factor = pow((l0/(2.0*mathematicalConstant::pi)),3.0);
  60
  61     y *= factor;
  62
  63     // and divide by the number of points in the box, to give the
  64     // energy density.
  65
  66     y /= scalar(K.size());
  67
  68     return kShellMeanEk;
  69 }
  70
  71
  72 // kShellMean : average over the points in a k-shell to evaluate the
  73 // radial part of the energy spectrum.
  74
  75 graph kShellMean
  76 (
  77     const complexVectorField& Ek,
  78     const Kmesh& K
  79 )
  80 {
  81     const label tnp = Ek.size();
  82     const label NoSubintervals = label
  83     (
  84         pow(scalar(tnp), 1.0/vector::dim)*pow(1.0/vector::dim, 0.5) - 0.5
  85     );
  86
  87     scalarField k1D(NoSubintervals);
  88     scalarField Ek1D(NoSubintervals);
  89     scalarField EWeight(NoSubintervals);
  90
  91     scalar kmax = K.max()*pow(1.0/vector::dim,0.5);
  92     scalar delta_k = kmax/(NoSubintervals);
  93
  94     forAll(Ek1D, a)
  95     {
  96         k1D[a] = (a + 1)*delta_k;
  97         Ek1D[a] = 0.0;
  98         EWeight[a] = 0;
  99     }
 100
 101     forAll(K, l)
 102     {
 103         scalar kmag = mag(K[l]);
 104
 105         for (label a=0; a<NoSubintervals; a++)
 106         {
 107             if
 108             (
 109                 kmag <= ((a + 1)*delta_k + delta_k/2.0)
 110              && kmag > ((a + 1)*delta_k - delta_k/2.0)
 111             )
 112             {
 113                 scalar dist = delta_k/2.0 - mag((a + 1)*delta_k - kmag);
 114
 115                 Ek1D[a] += dist*
 116                 magSqr
 117                 (
 118                     vector
 119                     (
 120                         mag(Ek[l].x()),
 121                         mag(Ek[l].y()),
 122                         mag(Ek[l].z())
 123                     )
 124                  );
 125
 126                 EWeight[a] += dist;
 127             }
 128         }
 129     }
 130
 131     for (label a=0; a<NoSubintervals; a++)
 132     {
 133         if (EWeight[a] > 0)
 134         {
 135             Ek1D[a] /= EWeight[a];
 136         }
 137     }
 138
 139     return graph("E(k)", "k", "E(k)", Ek1D, k1D);
 140 }
 141
 142
 143 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 144
 145 } // End namespace Foam
 146
 147 // ************************************************************************* //
 148