potential.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 "potential.H"
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
31 {
32  DynamicList<word> siteIdList;
33  DynamicList<word> pairPotentialSiteIdList;
34 
35  forAll(idList_, i)
36  {
37  const word& id(idList_[i]);
38 
39  if (!moleculePropertiesDict.found(id))
40  {
42  << id << " molecule subDict not found"
43  << nl << abort(FatalError);
44  }
45 
46  const dictionary& molDict(moleculePropertiesDict.subDict(id));
47 
48  List<word> siteIdNames = molDict.lookup("siteIds");
49 
50  forAll(siteIdNames, sI)
51  {
52  const word& siteId = siteIdNames[sI];
53 
54  if (findIndex(siteIdList, siteId) == -1)
55  {
56  siteIdList.append(siteId);
57  }
58  }
59 
60  List<word> pairPotSiteIds = molDict.lookup("pairPotentialSiteIds");
61 
62  forAll(pairPotSiteIds, sI)
63  {
64  const word& siteId = pairPotSiteIds[sI];
65 
66  if (findIndex(siteIdNames, siteId) == -1)
67  {
69  << siteId << " in pairPotentialSiteIds is not in siteIds: "
70  << siteIdNames << nl << abort(FatalError);
71  }
72 
73  if (findIndex(pairPotentialSiteIdList, siteId) == -1)
74  {
75  pairPotentialSiteIdList.append(siteId);
76  }
77  }
78  }
79 
80  nPairPotIds_ = pairPotentialSiteIdList.size();
81 
82  forAll(siteIdList, aSIN)
83  {
84  const word& siteId = siteIdList[aSIN];
85 
86  if (findIndex(pairPotentialSiteIdList, siteId) == -1)
87  {
88  pairPotentialSiteIdList.append(siteId);
89  }
90  }
91 
92  siteIdList_.transfer(pairPotentialSiteIdList.shrink());
93 }
94 
95 
96 void Foam::potential::potential::readPotentialDict()
97 {
98  Info<< nl << "Reading potential dictionary:" << endl;
99 
100  IOdictionary idListDict
101  (
102  IOobject
103  (
104  "idList",
105  mesh_.time().constant(),
106  mesh_,
109  )
110  );
111 
112  idList_ = List<word>(idListDict.lookup("idList"));
113 
114  setSiteIdList
115  (
116  IOdictionary
117  (
118  IOobject
119  (
120  "moleculeProperties",
121  mesh_.time().constant(),
122  mesh_,
125  false
126  )
127  )
128  );
129 
130  List<word> pairPotentialSiteIdList
131  (
132  SubList<word>(siteIdList_, nPairPotIds_)
133  );
134 
135  Info<< nl << "Unique site ids found: " << siteIdList_
136  << nl << "Site Ids requiring a pair potential: "
137  << pairPotentialSiteIdList
138  << endl;
139 
140  List<word> tetherSiteIdList(0);
141 
142  if (idListDict.found("tetherSiteIdList"))
143  {
144  tetherSiteIdList = List<word>(idListDict.lookup("tetherSiteIdList"));
145  }
146 
147  IOdictionary potentialDict
148  (
149  IOobject
150  (
151  "potentialDict",
152  mesh_.time().system(),
153  mesh_,
156  )
157  );
158 
159  potentialEnergyLimit_ =
160  potentialDict.lookup<scalar>("potentialEnergyLimit");
161 
162  if (potentialDict.found("removalOrder"))
163  {
164  List<word> remOrd = potentialDict.lookup("removalOrder");
165 
166  removalOrder_.setSize(remOrd.size());
167 
168  forAll(removalOrder_, rO)
169  {
170  removalOrder_[rO] = findIndex(idList_, remOrd[rO]);
171 
172  if (removalOrder_[rO] == -1)
173  {
175  << "removalOrder entry: " << remOrd[rO]
176  << " not found in idList."
177  << nl << abort(FatalError);
178  }
179  }
180  }
181 
182  // *************************************************************************
183  // Pair potentials
184 
185  if (!potentialDict.found("pair"))
186  {
188  << "pair potential specification subDict not found"
189  << abort(FatalError);
190  }
191 
192  const dictionary& pairDict = potentialDict.subDict("pair");
193 
194  pairPotentials_.buildPotentials
195  (
196  pairPotentialSiteIdList,
197  pairDict,
198  mesh_
199  );
200 
201  // *************************************************************************
202  // Tether potentials
203 
204  if (tetherSiteIdList.size())
205  {
206  if (!potentialDict.found("tether"))
207  {
209  << "tether potential specification subDict not found"
210  << abort(FatalError);
211  }
212 
213  const dictionary& tetherDict = potentialDict.subDict("tether");
214 
215  tetherPotentials_.buildPotentials
216  (
217  siteIdList_,
218  tetherDict,
219  tetherSiteIdList
220  );
221  }
222 
223  // *************************************************************************
224  // External Forces
225 
226  gravity_ = Zero;
227 
228  if (potentialDict.found("external"))
229  {
230  Info<< nl << "Reading external forces:" << endl;
231 
232  const dictionary& externalDict = potentialDict.subDict("external");
233 
234  // gravity
235  externalDict.readIfPresent("gravity", gravity_);
236  }
237 
238  Info<< nl << tab << "gravity = " << gravity_ << endl;
239 }
240 
241 
242 void Foam::potential::potential::readMdInitialiseDict
243 (
244  const IOdictionary& mdInitialiseDict,
245  IOdictionary& idListDict
246 )
247 {
248  IOdictionary moleculePropertiesDict
249  (
250  IOobject
251  (
252  "moleculeProperties",
253  mesh_.time().constant(),
254  mesh_,
257  false
258  )
259  );
260 
261  DynamicList<word> idList;
262 
263  DynamicList<word> tetherSiteIdList;
264 
265  forAll(mdInitialiseDict.toc(), zone)
266  {
267  const dictionary& zoneDict = mdInitialiseDict.subDict
268  (
269  mdInitialiseDict.toc()[zone]
270  );
271 
272  List<word> latticeIds
273  (
274  zoneDict.lookup("latticeIds")
275  );
276 
277  forAll(latticeIds, i)
278  {
279  const word& id = latticeIds[i];
280 
281  if (!moleculePropertiesDict.found(id))
282  {
284  << "Molecule type " << id
285  << " not found in moleculeProperties dictionary." << nl
286  << abort(FatalError);
287  }
288 
289  if (findIndex(idList,id) == -1)
290  {
291  idList.append(id);
292  }
293  }
294 
295  List<word> tetherSiteIds
296  (
297  zoneDict.lookup("tetherSiteIds")
298  );
299 
300  forAll(tetherSiteIds, t)
301  {
302  const word& tetherSiteId = tetherSiteIds[t];
303 
304  bool idFound = false;
305 
306  forAll(latticeIds, i)
307  {
308  if (idFound)
309  {
310  break;
311  }
312 
313  const word& id = latticeIds[i];
314 
315  List<word> siteIds
316  (
317  moleculePropertiesDict.subDict(id).lookup("siteIds")
318  );
319 
320  if (findIndex(siteIds, tetherSiteId) != -1)
321  {
322  idFound = true;
323  }
324  }
325 
326  if (idFound)
327  {
328  tetherSiteIdList.append(tetherSiteId);
329  }
330  else
331  {
333  << " not found as a site of any molecule in zone." << nl
334  << abort(FatalError);
335  }
336  }
337  }
338 
339  idList_.transfer(idList);
340 
341  tetherSiteIdList.shrink();
342 
343  idListDict.add("idList", idList_);
344 
345  idListDict.add("tetherSiteIdList", tetherSiteIdList);
346 
347  setSiteIdList(moleculePropertiesDict);
348 }
349 
350 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
351 
353 :
354  mesh_(mesh)
355 {
356  readPotentialDict();
357 }
358 
359 
361 (
362  const polyMesh& mesh,
363  const IOdictionary& mdInitialiseDict,
364  IOdictionary& idListDict
365 )
366 :
367  mesh_(mesh)
368 {
369  readMdInitialiseDict(mdInitialiseDict, idListDict);
370 }
371 
372 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
373 
375 {}
376 
377 
378 // ************************************************************************* //
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:663
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
static const char tab
Definition: Ostream.H:259
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
~potential()
Destructor.
Definition: potential.C:374
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
const word & constant() const
Return constant name.
Definition: TimePaths.H:123
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const word & system() const
Return system name.
Definition: TimePaths.H:113
const List< word > & idList() const
Definition: potentialI.H:34
static const char nl
Definition: Ostream.H:260
const Time & time() const
Return time.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void buildPotentials(const List< word > &siteIdList, const dictionary &tetherPotentialDict, const List< word > &tetherSiteIdList)
potential(const polyMesh &mesh)
Construct from mesh reference.
Definition: potential.C:352
void setSize(const label)
Reset size of List.
Definition: List.C:281
void buildPotentials(const List< word > &idList, const dictionary &pairPotentialDict, const polyMesh &mesh)
messageStream Info
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864
const List< word > & siteIdList() const
Definition: potentialI.H:40