src/randomProcesses/fft/kShellIntegration.C

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