src/randomProcesses/fft/kShellIntegration.C

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