pairPotential.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
39 
40 void Foam::pairPotential::scaleEnergy(scalar& e, const scalar r) const
41 {
42  if (!esfPtr_)
43  {
45  (
46  name_, pairPotentialProperties_, *this
47  ).ptr();
48  }
49 
50  esfPtr_->scaleEnergy(e, r);
51 }
52 
53 
54 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
55 
57 (
58  const word& name,
59  const dictionary& pairPotentialProperties
60 )
61 :
62  name_(name),
63  pairPotentialProperties_(pairPotentialProperties),
64  rCut_(pairPotentialProperties_.template lookup<scalar>("rCut")),
65  rCutSqr_(rCut_*rCut_),
66  rMin_(pairPotentialProperties_.template lookup<scalar>("rMin")),
67  dr_(pairPotentialProperties_.template lookup<scalar>("dr")),
68  forceLookup_(0),
69  energyLookup_(0),
70  esfPtr_(nullptr),
71  writeTables_(Switch(pairPotentialProperties_.lookup("writeTables")))
72 {}
73 
74 
75 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
76 
78 {
79  label N = label((rCut_ - rMin_)/dr_) + 1;
80 
81  forceLookup_.setSize(N);
82 
83  energyLookup_.setSize(N);
84 
85  forAll(forceLookup_, k)
86  {
87  energyLookup_[k] = scaledEnergy(k*dr_ + rMin_);
88 
89  forceLookup_[k] = -energyDerivative((k*dr_ + rMin_), true);
90  }
91 }
92 
93 
94 Foam::scalar Foam::pairPotential::force(const scalar r) const
95 {
96  scalar k_rIJ = (r - rMin_)/dr_;
97 
98  label k = label(k_rIJ);
99 
100  if (k < 0)
101  {
103  << "r less than rMin in pair potential " << name_ << nl
104  << abort(FatalError);
105  }
106 
107  scalar f =
108  (k_rIJ - k)*forceLookup_[k+1]
109  + (k + 1 - k_rIJ)*forceLookup_[k];
110 
111  return f;
112 }
113 
114 
117 {
118  List<Pair<scalar>> forceTab(forceLookup_.size());
119 
120  forAll(forceLookup_,k)
121  {
122  forceTab[k].first() = rMin_ + k*dr_;
123 
124  forceTab[k].second() = forceLookup_[k];
125  }
126 
127  return forceTab;
128 }
129 
130 
131 Foam::scalar Foam::pairPotential::energy(const scalar r) const
132 {
133  scalar k_rIJ = (r - rMin_)/dr_;
134 
135  label k = label(k_rIJ);
136 
137  if (k < 0)
138  {
140  << "r less than rMin in pair potential " << name_ << nl
141  << abort(FatalError);
142  }
143 
144  scalar e =
145  (k_rIJ - k)*energyLookup_[k+1]
146  + (k + 1 - k_rIJ)*energyLookup_[k];
147 
148  return e;
149 }
150 
151 
154 {
155  List<Pair<scalar>> energyTab(energyLookup_.size());
156 
157  forAll(energyLookup_,k)
158  {
159  energyTab[k].first() = rMin_ + k*dr_;
160 
161  energyTab[k].second() = energyLookup_[k];
162  }
163 
164  return energyTab;
165 }
166 
167 
168 Foam::scalar Foam::pairPotential::scaledEnergy(const scalar r) const
169 {
170  scalar e = unscaledEnergy(r);
171 
172  scaleEnergy(e, r);
173 
174  return e;
175 }
176 
177 
179 (
180  const scalar r,
181  const bool scaledEnergyDerivative
182 ) const
183 {
184  // Local quadratic fit to energy: E = a0 + a1*r + a2*r^2
185  // Differentiate to give f = -dE/dr = -a1 - 2*a2*r
186 
187  scalar ra = r - dr_;
188  scalar rf = r;
189  scalar rb = r + dr_;
190 
191  scalar Ea, Ef, Eb;
192 
193  if (scaledEnergyDerivative)
194  {
195  Ea = scaledEnergy(ra);
196  Ef = scaledEnergy(rf);
197  Eb = scaledEnergy(rb);
198  }
199  else
200  {
201  Ea = unscaledEnergy(ra);
202  Ef = unscaledEnergy(rf);
203  Eb = unscaledEnergy(rb);
204  }
205 
206  scalar denominator = (ra - rf)*(ra - rb)*(rf - rb);
207 
208  scalar a1 =
209  (
210  rb*rb*(Ea - Ef) + ra*ra*(Ef - Eb) + rf*rf*(Eb - Ea)
211  ) / denominator;
212 
213  scalar a2 =
214  (
215  rb*(Ef - Ea) + rf*(Ea - Eb) + ra*(Eb - Ef)
216  ) / denominator;
217 
218  return a1 + 2.0*a2*r;
219 }
220 
221 
222 bool Foam::pairPotential::read(const dictionary& pairPotentialProperties)
223 {
224  pairPotentialProperties_ = pairPotentialProperties;
225 
226  return true;
227 }
228 
229 
230 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
static autoPtr< energyScalingFunction > New(const word &name, const dictionary &energyScalingFunctionProperties, const pairPotential &pairPot)
Return a reference to the selected viscosity model.
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
virtual bool read(const dictionary &pairPotentialProperties)=0
Read pairPotential dictionary.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
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
scalar energyDerivative(const scalar r, const bool scaledEnergyDerivative=true) 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/any.
Definition: Switch.H:60
label k
Boltzmann constant.
List< Pair< scalar > > forceTable() const
scalar energy(const scalar r) const
pairPotential(const word &name, const dictionary &pairPotentialProperties)
Construct from components.
Definition: pairPotential.C:57
A class for handling words, derived from string.
Definition: word.H:59
scalar force(const scalar r) const
Definition: pairPotential.C:94
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static const char nl
Definition: Ostream.H:260
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
List< Pair< scalar > > energyTable() const
void scaleEnergy(scalar &e, const scalar r) const
Definition: pairPotential.C:40
label N
Definition: createFields.H:22
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
scalar scaledEnergy(const scalar r) const
scalar Ea(const scalar p, const scalar T) const
Definition: HtoEthermo.H:20
Namespace for OpenFOAM.