pairPotential.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 "pairPotential.H"
27 #include "energyScalingFunction.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(pairPotential, 0);
34  defineRunTimeSelectionTable(pairPotential, dictionary);
35 }
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 
41 void Foam::pairPotential::scaleEnergy(scalar& e, const scalar r) const
42 {
43  if (!esfPtr_)
44  {
46  (
47  name_, pairPotentialProperties_, *this
48  ).ptr();
49  }
50 
51  esfPtr_->scaleEnergy(e, r);
52 }
53 
54 
55 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
56 
58 (
59  const word& name,
60  const dictionary& pairPotentialProperties
61 )
62 :
63  name_(name),
64  pairPotentialProperties_(pairPotentialProperties),
65  rCut_(readScalar(pairPotentialProperties_.lookup("rCut"))),
66  rCutSqr_(rCut_*rCut_),
67  rMin_(readScalar(pairPotentialProperties_.lookup("rMin"))),
68  dr_(readScalar(pairPotentialProperties_.lookup("dr"))),
69  forceLookup_(0),
70  energyLookup_(0),
71  esfPtr_(NULL),
72  writeTables_(Switch(pairPotentialProperties_.lookup("writeTables")))
73 {}
74 
75 
76 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
77 
79 {
80  label N = label((rCut_ - rMin_)/dr_) + 1;
81 
82  forceLookup_.setSize(N);
83 
84  energyLookup_.setSize(N);
85 
86  forAll(forceLookup_, k)
87  {
88  energyLookup_[k] = scaledEnergy(k*dr_ + rMin_);
89 
90  forceLookup_[k] = -energyDerivative((k*dr_ + rMin_), true);
91  }
92 }
93 
94 
95 Foam::scalar Foam::pairPotential::force(const scalar r) const
96 {
97  scalar k_rIJ = (r - rMin_)/dr_;
98 
99  label k = label(k_rIJ);
100 
101  if (k < 0)
102  {
103  FatalErrorIn("pairPotential::force(const scalar) const")
104  << "r less than rMin in pair potential " << name_ << nl
105  << abort(FatalError);
106  }
107 
108  scalar f =
109  (k_rIJ - k)*forceLookup_[k+1]
110  + (k + 1 - k_rIJ)*forceLookup_[k];
111 
112  return f;
113 }
114 
115 
118 {
119  List<Pair<scalar> > forceTab(forceLookup_.size());
120 
121  forAll(forceLookup_,k)
122  {
123  forceTab[k].first() = rMin_ + k*dr_;
124 
125  forceTab[k].second() = forceLookup_[k];
126  }
127 
128  return forceTab;
129 }
130 
131 
132 Foam::scalar Foam::pairPotential::energy(const scalar r) const
133 {
134  scalar k_rIJ = (r - rMin_)/dr_;
135 
136  label k = label(k_rIJ);
137 
138  if (k < 0)
139  {
140  FatalErrorIn("pairPotential::energy(const scalar) const")
141  << "r less than rMin in pair potential " << name_ << nl
142  << abort(FatalError);
143  }
144 
145  scalar e =
146  (k_rIJ - k)*energyLookup_[k+1]
147  + (k + 1 - k_rIJ)*energyLookup_[k];
148 
149  return e;
150 }
151 
152 
155 {
156  List<Pair<scalar> > energyTab(energyLookup_.size());
157 
158  forAll(energyLookup_,k)
159  {
160  energyTab[k].first() = rMin_ + k*dr_;
161 
162  energyTab[k].second() = energyLookup_[k];
163  }
164 
165  return energyTab;
166 }
167 
168 
169 Foam::scalar Foam::pairPotential::scaledEnergy(const scalar r) const
170 {
171  scalar e = unscaledEnergy(r);
172 
173  scaleEnergy(e, r);
174 
175  return e;
176 }
177 
178 
180 (
181  const scalar r,
182  const bool scaledEnergyDerivative
183 ) const
184 {
185  // Local quadratic fit to energy: E = a0 + a1*r + a2*r^2
186  // Differentiate to give f = -dE/dr = -a1 - 2*a2*r
187 
188  scalar ra = r - dr_;
189  scalar rf = r;
190  scalar rb = r + dr_;
191 
192  scalar Ea, Ef, Eb;
193 
194  if (scaledEnergyDerivative)
195  {
196  Ea = scaledEnergy(ra);
197  Ef = scaledEnergy(rf);
198  Eb = scaledEnergy(rb);
199  }
200  else
201  {
202  Ea = unscaledEnergy(ra);
203  Ef = unscaledEnergy(rf);
204  Eb = unscaledEnergy(rb);
205  }
206 
207  scalar denominator = (ra - rf)*(ra - rb)*(rf - rb);
208 
209  scalar a1 =
210  (
211  rb*rb*(Ea - Ef) + ra*ra*(Ef - Eb) + rf*rf*(Eb - Ea)
212  ) / denominator;
213 
214  scalar a2 =
215  (
216  rb*(Ef - Ea) + rf*(Ea - Eb) + ra*(Eb - Ef)
217  ) / denominator;
218 
219  return a1 + 2.0*a2*r;
220 }
221 
222 
223 bool Foam::pairPotential::read(const dictionary& pairPotentialProperties)
224 {
225  pairPotentialProperties_ = pairPotentialProperties;
226 
227  return true;
228 }
229 
230 
231 // ************************************************************************* //
void scaleEnergy(scalar &e, const scalar r) const
Definition: pairPotential.C:41
scalar scaledEnergy(const scalar r) const
virtual bool read(const dictionary &pairPotentialProperties)=0
Read pairPotential dictionary.
labelList f(nPoints)
A class for handling words, derived from string.
Definition: word.H:59
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
List< Pair< scalar > > forceTable() const
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
pairPotential(const pairPotential &)
Disallow copy construct.
scalar energyDerivative(const scalar r, const bool scaledEnergyDerivative=true) const
Namespace for OpenFOAM.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
static const char nl
Definition: Ostream.H:260
const double e
Elementary charge.
Definition: doubleFloat.H:78
scalar force(const scalar r) const
Definition: pairPotential.C:95
scalar energy(const scalar r) const
label N
Definition: createFields.H:22
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
label k
Boltzmann constant.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< Pair< scalar > > energyTable() const
error FatalError
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
static autoPtr< energyScalingFunction > New(const word &name, const dictionary &energyScalingFunctionProperties, const pairPotential &pairPot)
Return a reference to the selected viscosity model.
defineTypeNameAndDebug(combustionModel, 0)