kShellIntegration.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
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 
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 
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 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Class to create, store and output qgraph files.
Definition: graph.H:58
scalar y
graph kShellMean(const complexVectorField &Ek, const Kmesh &K)
iterator begin()
Iterator set to the beginning of the HashTable.
Definition: HashTableI.H:419
Calculate the wavenumber vector field corresponding to the space vector field of a finite volume mesh...
Definition: Kmesh.H:49
const vector & sizeOfBox() const
Definition: Kmesh.H:89
scalar max() const
Definition: Kmesh.H:99
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Integrate a multi-dimensional complexVectorField in k-shells to create the 1D.
graph kShellIntegration(const complexVectorField &Ek, const Kmesh &K)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const labelList & nn() const
Definition: Kmesh.H:94
dimensioned< scalar > mag(const dimensioned< Type > &)
const scalarField & x() const
Definition: graph.H:165