pairPotentialList.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2019-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "pairPotentialList.H"
30 #include "OFstream.H"
31 #include "Time.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 void Foam::pairPotentialList::readPairPotentialDict
36 (
37  const List<word>& idList,
38  const dictionary& pairPotentialDict,
39  const polyMesh& mesh
40 )
41 {
42  Info<< nl << "Building pair potentials." << endl;
43 
44  rCutMax_ = 0.0;
45 
46  for (label a = 0; a < nIds_; ++a)
47  {
48  word idA = idList[a];
49 
50  for (label b = a; b < nIds_; ++b)
51  {
52  word idB = idList[b];
53 
54  word pairPotentialName;
55 
56  if (a == b)
57  {
58  if (pairPotentialDict.found(idA + "-" + idB))
59  {
60  pairPotentialName = idA + "-" + idB;
61  }
62  else
63  {
65  << "Pair pairPotential specification subDict "
66  << idA << "-" << idB << " not found"
67  << nl << abort(FatalError);
68  }
69  }
70  else
71  {
72  if (pairPotentialDict.found(idA + "-" + idB))
73  {
74  pairPotentialName = idA + "-" + idB;
75  }
76 
77  else if (pairPotentialDict.found(idB + "-" + idA))
78  {
79  pairPotentialName = idB + "-" + idA;
80  }
81 
82  else
83  {
85  << "Pair pairPotential specification subDict "
86  << idA << "-" << idB << " or "
87  << idB << "-" << idA << " not found"
88  << nl << abort(FatalError);
89  }
90 
91  if
92  (
93  pairPotentialDict.found(idA+"-"+idB)
94  && pairPotentialDict.found(idB+"-"+idA)
95  )
96  {
98  << "Pair pairPotential specification subDict "
99  << idA << "-" << idB << " and "
100  << idB << "-" << idA << " found multiple definition"
101  << nl << abort(FatalError);
102  }
103  }
104 
105  (*this).set
106  (
107  pairPotentialIndex(a, b),
109  (
110  pairPotentialName,
111  pairPotentialDict.subDict(pairPotentialName)
112  )
113  );
114 
115  if ((*this)[pairPotentialIndex(a, b)].rCut() > rCutMax_)
116  {
117  rCutMax_ = (*this)[pairPotentialIndex(a, b)].rCut();
118  }
119 
120  if ((*this)[pairPotentialIndex(a, b)].writeTables())
121  {
122  fileHandler().mkDir(mesh.time().path());
123  autoPtr<OSstream> ppTabFile
124  (
125  fileHandler().NewOFstream
126  (
127  mesh.time().path()/pairPotentialName
128  )
129  );
130 
131  if
132  (
133  !(*this)[pairPotentialIndex(a, b)].writeEnergyAndForceTables
134  (
135  ppTabFile()
136  )
137  )
138  {
140  << "Failed writing to "
141  << ppTabFile().name() << nl
142  << abort(FatalError);
143  }
144  }
145  }
146  }
147 
148  if (!pairPotentialDict.found("electrostatic"))
149  {
151  << "Pair pairPotential specification subDict electrostatic"
152  << nl << abort(FatalError);
153  }
154 
155  electrostaticPotential_ = pairPotential::New
156  (
157  "electrostatic",
158  pairPotentialDict.subDict("electrostatic")
159  );
160 
161  if (electrostaticPotential_->rCut() > rCutMax_)
162  {
163  rCutMax_ = electrostaticPotential_->rCut();
164  }
165 
166  if (electrostaticPotential_->writeTables())
167  {
168  fileHandler().mkDir(mesh.time().path());
169  autoPtr<OSstream> ppTabFile
170  (
171  fileHandler().NewOFstream
172  (
173  mesh.time().path()/"electrostatic"
174  )
175  );
176 
177  if (!electrostaticPotential_->writeEnergyAndForceTables(ppTabFile()))
178  {
180  << "Failed writing to "
181  << ppTabFile().name() << nl
182  << abort(FatalError);
183  }
184  }
185 
186  rCutMaxSqr_ = rCutMax_*rCutMax_;
187 }
188 
189 
190 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
191 
193 :
195 {}
196 
197 
199 (
200  const List<word>& idList,
201  const dictionary& pairPotentialDict,
202  const polyMesh& mesh
203 )
204 :
206 {
207  buildPotentials(idList, pairPotentialDict, mesh);
208 }
209 
210 
211 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
212 
214 {}
215 
216 
217 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
218 
220 (
221  const List<word>& idList,
222  const dictionary& pairPotentialDict,
223  const polyMesh& mesh
224 )
225 {
226  setSize(((idList.size()*(idList.size() + 1))/2));
227 
228  nIds_ = idList.size();
229 
230  readPairPotentialDict(idList, pairPotentialDict, mesh);
231 }
232 
233 
235 (
236  const label a,
237  const label b
238 ) const
239 {
240  return (*this)[pairPotentialIndex(a, b)];
241 }
242 
243 
244 bool Foam::pairPotentialList::rCutMaxSqr(const scalar rIJMagSqr) const
245 {
246  if (rIJMagSqr < rCutMaxSqr_)
247  {
248  return true;
249  }
250 
251  return false;
252 }
253 
254 
256 (
257  const label a,
258  const label b,
259  const scalar rIJMagSqr
260 ) const
261 {
262  if (rIJMagSqr < rCutSqr(a, b))
263  {
264  return true;
265  }
266 
267  return false;
268 }
269 
270 
272 (
273  const label a,
274  const label b
275 ) const
276 {
277  return (*this)[pairPotentialIndex(a, b)].rMin();
278 }
279 
280 
281 Foam::scalar Foam::pairPotentialList::dr
282 (
283  const label a,
284  const label b
285 ) const
286 {
287  return (*this)[pairPotentialIndex(a, b)].dr();
288 }
289 
290 
292 (
293  const label a,
294  const label b
295 ) const
296 {
297  return (*this)[pairPotentialIndex(a, b)].rCutSqr();
298 }
299 
300 
302 (
303  const label a,
304  const label b
305 ) const
306 {
307  return (*this)[pairPotentialIndex(a, b)].rCut();
308 }
309 
310 
312 (
313  const label a,
314  const label b,
315  const scalar rIJMag
316 ) const
317 {
318  scalar f = (*this)[pairPotentialIndex(a, b)].force(rIJMag);
319 
320  return f;
321 }
322 
323 
325 (
326  const label a,
327  const label b,
328  const scalar rIJMag
329 ) const
330 {
331  scalar e = (*this)[pairPotentialIndex(a, b)].energy(rIJMag);
332 
333  return e;
334 }
335 
336 
337 // ************************************************************************* //
Foam::pairPotential::New
static autoPtr< pairPotential > New(const word &name, const dictionary &pairPotentialProperties)
Return a reference to the selected viscosity model.
Definition: pairPotentialNew.C:35
setSize
points setSize(newPointi)
Foam::pairPotentialList::rCut
scalar rCut(const label a, const label b) const
Definition: pairPotentialList.C:302
Foam::fileHandler
const fileOperation & fileHandler()
Get current file handler.
Definition: fileOperation.C:1485
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
Foam::pairPotentialList::rCutMaxSqr
scalar rCutMaxSqr() const
Definition: pairPotentialListI.H:68
OFstream.H
Foam::pairPotentialList::rCutSqr
bool rCutSqr(const label a, const label b, const scalar rIJMagSqr) const
Definition: pairPotentialList.C:256
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::pairPotentialList::energy
scalar energy(const label a, const label b, const scalar rIJMag) const
Definition: pairPotentialList.C:325
Foam::pairPotentialList::~pairPotentialList
~pairPotentialList()
Destructor.
Definition: pairPotentialList.C:213
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::pairPotentialList::rMin
scalar rMin(const label a, const label b) const
Definition: pairPotentialList.C:272
Foam::pairPotentialList::dr
scalar dr(const label a, const label b) const
Definition: pairPotentialList.C:282
Foam::PtrList
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: List.H:59
pairPotentialList.H
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::pairPotential
Definition: pairPotential.H:60
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::pairPotentialList::force
scalar force(const label a, const label b, const scalar rIJMag) const
Definition: pairPotentialList.C:312
Time.H
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::pairPotentialList::pairPotentialFunction
const pairPotential & pairPotentialFunction(const label a, const label b) const
Definition: pairPotentialList.C:235
Foam::nl
constexpr char nl
Definition: Ostream.H:404
f
labelList f(nPoints)
Foam::List< word >
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::fileOperation::mkDir
virtual bool mkDir(const fileName &, mode_t=0777) const =0
Make directory.
Foam::pairPotentialList::buildPotentials
void buildPotentials(const List< word > &idList, const dictionary &pairPotentialDict, const polyMesh &mesh)
Definition: pairPotentialList.C:220
Foam::pairPotentialList::pairPotentialList
pairPotentialList()
Definition: pairPotentialList.C:192