src/randomProcesses/processes/UOprocess/UOprocess.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "error.H"
  27
  28 #include "UOprocess.H"
  29 #include "Kmesh.H"
  30 #include "dictionary.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 complexVector UOprocess::WeinerProcess()
  40 {
  41     return RootDeltaT*complexVector
  42     (
  43         complex(GaussGen.GaussNormal(), GaussGen.GaussNormal()),
  44         complex(GaussGen.GaussNormal(), GaussGen.GaussNormal()),
  45         complex(GaussGen.GaussNormal(), GaussGen.GaussNormal())
  46     );
  47 }
  48
  49
  50 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  51
  52 // from components
  53 UOprocess::UOprocess
  54 (
  55     const Kmesh& kmesh,
  56     const scalar deltaT,
  57     const dictionary& UOdict
  58 )
  59 :
  60     GaussGen(label(0)),
  61     Mesh(kmesh),
  62     DeltaT(deltaT),
  63     RootDeltaT(sqrt(DeltaT)),
  64     UOfield(Mesh.size()),
  65
  66     Alpha(readScalar(UOdict.lookup("UOalpha"))),
  67     Sigma(readScalar(UOdict.lookup("UOsigma"))),
  68     Kupper(readScalar(UOdict.lookup("UOKupper"))),
  69     Klower(readScalar(UOdict.lookup("UOKlower"))),
  70     Scale((Kupper - Klower)*pow(scalar(Mesh.size()), 1.0/vector::dim))
  71 {
  72     const vectorField& K = Mesh;
  73
  74     scalar sqrKupper = sqr(Kupper);
  75     scalar sqrKlower = sqr(Klower) + SMALL;
  76     scalar sqrK;
  77
  78     forAll(UOfield, i)
  79     {
  80         if ((sqrK = magSqr(K[i])) < sqrKupper && sqrK > sqrKlower)
  81         {
  82             UOfield[i] = Scale*Sigma*WeinerProcess();
  83         }
  84         else
  85         {
  86             UOfield[i] = complexVector
  87             (
  88                 complex(0, 0),
  89                 complex(0, 0),
  90                 complex(0, 0)
  91             );
  92         }
  93     }
  94 }
  95
  96
  97 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  98
  99 const complexVectorField& UOprocess::newField()
 100 {
 101     const vectorField& K = Mesh;
 102
 103     label count = 0;
 104     scalar sqrKupper = sqr(Kupper);
 105     scalar sqrKlower = sqr(Klower) + SMALL;
 106     scalar sqrK;
 107
 108     forAll(UOfield, i)
 109     {
 110         if ((sqrK = magSqr(K[i])) < sqrKupper && sqrK > sqrKlower)
 111         {
 112             count++;
 113             UOfield[i] =
 114                 (1.0 - Alpha*DeltaT)*UOfield[i]
 115               + Scale*Sigma*WeinerProcess();
 116         }
 117     }
 118
 119     Info<< "    Number of forced K = " << count << nl;
 120
 121     return UOfield;
 122 }
 123
 124
 125 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 126
 127 } // End namespace Foam
 128
 129 // ************************************************************************* //
 130