src/randomProcesses/processes/UOprocess/UOprocess.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "error.H"
  28
  29 #include "UOprocess.H"
  30 #include "Kmesh.H"
  31 #include "dictionary.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 complexVector UOprocess::WeinerProcess()
  41 {
  42     return RootDeltaT*complexVector
  43     (
  44         complex(GaussGen.GaussNormal(), GaussGen.GaussNormal()),
  45         complex(GaussGen.GaussNormal(), GaussGen.GaussNormal()),
  46         complex(GaussGen.GaussNormal(), GaussGen.GaussNormal())
  47     );
  48 }
  49
  50
  51 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  52
  53 // from components
  54 UOprocess::UOprocess
  55 (
  56     const Kmesh& kmesh,
  57     const scalar deltaT,
  58     const dictionary& UOdict
  59 )
  60 :
  61     GaussGen(label(0)),
  62     Mesh(kmesh),
  63     DeltaT(deltaT),
  64     RootDeltaT(sqrt(DeltaT)),
  65     UOfield(Mesh.size()),
  66
  67     Alpha(readScalar(UOdict.lookup("UOalpha"))),
  68     Sigma(readScalar(UOdict.lookup("UOsigma"))),
  69     Kupper(readScalar(UOdict.lookup("UOKupper"))),
  70     Klower(readScalar(UOdict.lookup("UOKlower"))),
  71     Scale((Kupper - Klower)*pow(scalar(Mesh.size()), 1.0/vector::dim))
  72 {
  73     const vectorField& K = Mesh;
  74
  75     scalar sqrKupper = sqr(Kupper);
  76     scalar sqrKlower = sqr(Klower) + SMALL;
  77     scalar sqrK;
  78
  79     forAll(UOfield, i)
  80     {
  81         if ((sqrK = magSqr(K[i])) < sqrKupper && sqrK > sqrKlower)
  82         {
  83             UOfield[i] = Scale*Sigma*WeinerProcess();
  84         }
  85         else
  86         {
  87             UOfield[i] = complexVector
  88             (
  89                 complex(0, 0),
  90                 complex(0, 0),
  91                 complex(0, 0)
  92             );
  93         }
  94     }
  95 }
  96
  97
  98 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  99
 100 const complexVectorField& UOprocess::newField()
 101 {
 102     const vectorField& K = Mesh;
 103
 104     label count = 0;
 105     scalar sqrKupper = sqr(Kupper);
 106     scalar sqrKlower = sqr(Klower) + SMALL;
 107     scalar sqrK;
 108
 109     forAll(UOfield, i)
 110     {
 111         if ((sqrK = magSqr(K[i])) < sqrKupper && sqrK > sqrKlower)
 112         {
 113             count++;
 114             UOfield[i] =
 115                 (1.0 - Alpha*DeltaT)*UOfield[i]
 116               + Scale*Sigma*WeinerProcess();
 117         }
 118     }
 119
 120     Info<< "    Number of forced K = " << count << nl;
 121
 122     return UOfield;
 123 }
 124
 125
 126 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 127
 128 } // End namespace Foam
 129
 130 // ************************************************************************* //
 131