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