zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry
2011-07-22 14:12:10 +01:00
parent 5f9e3a04f0 6be1a5d417
commit cc272b762e
26 changed files with 4 additions and 8084 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1265
src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
View File

File diff suppressed because it is too large Load Diff

252
src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
View File

@ -1,252 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::MoleculeCloud
Description
SourceFiles
MoleculeCloudI.H
MoleculeCloud.C
\*---------------------------------------------------------------------------*/
#ifndef MoleculeCloud_H
#define MoleculeCloud_H
#include "Cloud.H"
#include "moleculeCloud.H"
#include "IOdictionary.H"
#include "potential.H"
#include "InteractionLists.H"
#include "labelVector.H"
#include "Random.H"
#include "fileName.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class MoleculeCloud Declaration
\*---------------------------------------------------------------------------*/
template<class MoleculeType>
class MoleculeCloud
:
public Cloud<MoleculeType>,
public moleculeCloud
{
private:
// Private data
//-
const polyMesh& mesh_;
//-
const potential& pot_;
//-
List<DynamicList<MoleculeType*> > cellOccupancy_;
//-
InteractionLists<MoleculeType> il_;
//-
List<typename MoleculeType::constantProperties> constPropList_;
//-
Random rndGen_;
// Private Member Functions
//-
void buildConstProps();
//-
void setSiteSizesAndPositions();
//- Determine which molecules are in which cells
void buildCellOccupancy();
//-
void calculatePairForce();
//-
inline void evaluatePair
(
MoleculeType& molI,
MoleculeType& molJ
);
//-
inline bool evaluatePotentialLimit
(
MoleculeType& molI,
MoleculeType& molJ
) const;
//-
void calculateTetherForce();
//-
void calculateExternalForce();
//-
void removeHighEnergyOverlaps();
//-
void initialiseMolecules(const dictionary& mdInitialiseDict);
//-
void createMolecule
(
const point& position,
label cell,
label tetFace,
label tetPt,
label id,
bool tethered,
scalar temperature,
const vector& bulkVelocity
);
//-
label nSites() const;
//- Disallow default bitwise copy construct
MoleculeCloud(const MoleculeCloud&);
//- Disallow default bitwise assignment
void operator=(const MoleculeCloud&);
public:
// Constructors
//- Construct given mesh and potential references
MoleculeCloud
(
const word& cloudName,
const polyMesh& mesh,
const potential& pot,
bool readFields = true
);
//- Construct given mesh, potential and mdInitialiseDict
MoleculeCloud
(
const word& cloudName,
const polyMesh& mesh,
const potential& pot,
const dictionary& mdInitialiseDict,
bool readFields = true
);
// Member Functions
//- Evolve the molecules (move, calculate forces, control state etc)
void evolve();
//-
void calculateForce();
//- Print cloud information
void info();
// Access
//-
inline const polyMesh& mesh() const;
//-
inline const potential& pot() const;
//-
inline const List<DynamicList<MoleculeType*> >&
cellOccupancy() const;
//-
inline const InteractionLists<MoleculeType>& il() const;
//-
inline const List<typename MoleculeType::constantProperties>
constProps() const;
//-
inline const typename MoleculeType::constantProperties&
constProps(label id) const;
//-
inline Random& rndGen();
//-
inline vector equipartitionLinearVelocity
(
scalar temperature,
scalar mass
);
//-
inline vector equipartitionAngularMomentum
(
scalar temperature,
const typename MoleculeType::constantProperties& cP
);
// Member Operators
//- Write molecule sites in XYZ format
void writeXYZ(const fileName& fName) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "MoleculeCloudI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "MoleculeCloud.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

416
src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
View File

@ -1,416 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class MoleculeType>
inline void Foam::MoleculeCloud<MoleculeType>::evaluatePair
(
MoleculeType& molI,
MoleculeType& molJ
)
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI.id();
label idJ = molJ.id();
const typename MoleculeType::constantProperties& constPropI
(
constProps(idI)
);
const typename MoleculeType::constantProperties& constPropJ
(
constProps(idJ)
);
forAll(constPropI.pairPotSites(), pI)
{
label sI = constPropI.pairPotSites()[pI];
label idsI = constPropI.sites()[sI].siteId();
forAll(constPropJ.pairPotSites(), pJ)
{
label sJ = constPropJ.pairPotSites()[pJ];
label idsJ = constPropJ.sites()[sJ].siteId();
vector rsIsJ =
molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ);
if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
{
scalar rsIsJMag = mag(rsIsJ);
vector fsIsJ =
(rsIsJ/rsIsJMag)
*pairPot.force(idsI, idsJ, rsIsJMag);
molI.siteForces()[sI] += fsIsJ;
molJ.siteForces()[sJ] += -fsIsJ;
scalar potentialEnergy
(
pairPot.energy(idsI, idsJ, rsIsJMag)
);
molI.potentialEnergy() += 0.5*potentialEnergy;
molJ.potentialEnergy() += 0.5*potentialEnergy;
vector rIJ = molI.position() - molJ.position();
tensor virialContribution =
(rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
molI.rf() += virialContribution;
molJ.rf() += virialContribution;
}
}
}
forAll(constPropI.electrostaticSites(), pI)
{
label sI = constPropI.electrostaticSites()[pI];
forAll(constPropJ.electrostaticSites(), pJ)
{
label sJ = constPropJ.electrostaticSites()[pJ];
vector rsIsJ =
molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ);
if (rsIsJMagSq <= electrostatic.rCutSqr())
{
scalar rsIsJMag = mag(rsIsJ);
scalar chargeI = constPropI.sites()[sI].siteCharge();
scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
vector fsIsJ =
(rsIsJ/rsIsJMag)
*chargeI*chargeJ*electrostatic.force(rsIsJMag);
molI.siteForces()[sI] += fsIsJ;
molJ.siteForces()[sJ] += -fsIsJ;
scalar potentialEnergy =
chargeI*chargeJ
*electrostatic.energy(rsIsJMag);
molI.potentialEnergy() += 0.5*potentialEnergy;
molJ.potentialEnergy() += 0.5*potentialEnergy;
vector rIJ = molI.position() - molJ.position();
tensor virialContribution =
(rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
molI.rf() += virialContribution;
molJ.rf() += virialContribution;
}
}
}
}
template<class MoleculeType>
inline bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
(
MoleculeType& molI,
MoleculeType& molJ
) const
{
const pairPotentialList& pairPot = pot_.pairPotentials();
const pairPotential& electrostatic = pairPot.electrostatic();
label idI = molI.id();
label idJ = molJ.id();
const typename MoleculeType::constantProperties& constPropI
(
constProps(idI)
);
const typename MoleculeType::constantProperties& constPropJ
(
constProps(idJ)
);
forAll(constPropI.pairPotSites(), pI)
{
label sI = constPropI.pairPotSites()[pI];
label idsI = constPropI.sites()[sI].siteId();
forAll(constPropJ.pairPotSites(), pJ)
{
label sJ = constPropJ.pairPotSites()[pJ];
label idsJ = constPropJ.sites()[sJ].siteId();
vector rsIsJ =
molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ);
if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
{
scalar rsIsJMag = mag(rsIsJ);
// Guard against pairPot.energy being evaluated
// if rIJMag < SMALL. A floating point exception will
// happen otherwise.
if (rsIsJMag < SMALL)
{
WarningIn
(
"MoleculeCloud<MoleculeType>::"
"removeHighEnergyOverlaps()"
)
<< "Molecule site pair closer than "
<< SMALL
<< ": mag separation = " << rsIsJMag
<< ". These may have been placed on top of each"
<< " other by a rounding error in mdInitialise in"
<< " parallel or a block filled with moleculess"
<< " twice. Removing one of the molecules."
<< endl;
return true;
}
// Guard against pairPot.energy being evaluated if rIJMag <
// rMin. A tabulation lookup error will occur otherwise.
if (rsIsJMag < pairPot.rMin(idsI, idsJ))
{
return true;
}
if
(
mag(pairPot.energy(idsI, idsJ, rsIsJMag))
> pot_.potentialEnergyLimit()
)
{
return true;
};
}
}
}
forAll(constPropI.electrostaticSites(), pI)
{
label sI = constPropI.electrostaticSites()[pI];
forAll(constPropJ.electrostaticSites(), pJ)
{
label sJ = constPropJ.electrostaticSites()[pJ];
vector rsIsJ =
molI.sitePositions()[sI] - molJ.sitePositions()[sJ];
scalar rsIsJMagSq = magSqr(rsIsJ);
if (pairPot.rCutMaxSqr(rsIsJMagSq))
{
scalar rsIsJMag = mag(rsIsJ);
// Guard against pairPot.energy being evaluated
// if rIJMag < SMALL. A floating point exception will
// happen otherwise.
if (rsIsJMag < SMALL)
{
WarningIn
(
"MoleculeCloud<MoleculeType>::"
"removeHighEnergyOverlaps()"
)
<< "Molecule site pair closer than "
<< SMALL
<< ": mag separation = " << rsIsJMag
<< ". These may have been placed on top of each"
<< " other by a rounding error in mdInitialise in"
<< " parallel or a block filled with molecules"
<< " twice. Removing one of the molecules."
<< endl;
return true;
}
if (rsIsJMag < electrostatic.rMin())
{
return true;
}
scalar chargeI = constPropI.sites()[sI].siteCharge();
scalar chargeJ = constPropJ.sites()[sJ].siteCharge();
if
(
mag(chargeI*chargeJ*electrostatic.energy(rsIsJMag))
> pot_.potentialEnergyLimit()
)
{
return true;
};
}
}
}
return false;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class MoleculeType>
inline const Foam::polyMesh& Foam::MoleculeCloud<MoleculeType>::mesh() const
{
return mesh_;
}
template<class MoleculeType>
inline const Foam::potential& Foam::MoleculeCloud<MoleculeType>::pot() const
{
return pot_;
}
template<class MoleculeType>
inline const Foam::List<Foam::DynamicList<MoleculeType*> >&
Foam::MoleculeCloud<MoleculeType>::cellOccupancy() const
{
return cellOccupancy_;
}
template<class MoleculeType>
inline const Foam::InteractionLists<MoleculeType>&
Foam::MoleculeCloud<MoleculeType>::il() const
{
return il_;
}
template<class MoleculeType>
inline const Foam::List<typename MoleculeType::constantProperties>
Foam::MoleculeCloud<MoleculeType>::constProps() const
{
return constPropList_;
}
template<class MoleculeType>
inline const typename MoleculeType::constantProperties&
Foam::MoleculeCloud<MoleculeType>::constProps(label id) const
{
return constPropList_[id];
}
template<class MoleculeType>
inline Foam::Random& Foam::MoleculeCloud<MoleculeType>::rndGen()
{
return rndGen_;
}
template<class MoleculeType>
inline Foam::vector
Foam::MoleculeCloud<MoleculeType>::equipartitionLinearVelocity
(
scalar temperature,
scalar mass
)
{
return sqrt(physicoChemical::k.value()*temperature/mass)*vector
(
rndGen_.GaussNormal(),
rndGen_.GaussNormal(),
rndGen_.GaussNormal()
);
}
template<class MoleculeType>
inline Foam::vector
Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
(
scalar temperature,
const typename MoleculeType::constantProperties& cP
)
{
scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
if (cP.linearMolecule())
{
return sqrtKbT*vector
(
0.0,
sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
);
}
else
{
return sqrtKbT*vector
(
sqrt(cP.momentOfInertia().xx())*rndGen_.GaussNormal(),
sqrt(cP.momentOfInertia().yy())*rndGen_.GaussNormal(),
sqrt(cP.momentOfInertia().zz())*rndGen_.GaussNormal()
);
}
}
// ************************************************************************* //

48
src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.C
View File

@ -1,48 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "moleculeCloud.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(moleculeCloud, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::moleculeCloud::moleculeCloud()
{}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * //
Foam::moleculeCloud::~moleculeCloud()
{}
// ************************************************************************* //

83
src/lagrangian/molecularDynamics/molecule/clouds/baseClasses/moleculeCloud/moleculeCloud.H
View File

@ -1,83 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::moleculeCloud
Description
Virtual abstract base class for templated moleculeCloud
SourceFiles
moleculeCloud.C
\*---------------------------------------------------------------------------*/
#ifndef moleculeCloud_H
#define moleculeCloud_H
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class moleculeCloud Declaration
\*---------------------------------------------------------------------------*/
class moleculeCloud
{
// Private Member Functions
//- Disallow default bitwise copy construct
moleculeCloud(const moleculeCloud&);
//- Disallow default bitwise assignment
void operator=(const moleculeCloud&);
public:
//- Runtime type information
TypeName("moleculeCloud");
// Constructors
//- Null constructor
moleculeCloud();
//- Destructor
virtual ~moleculeCloud();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

52
src/lagrangian/molecularDynamics/molecule/clouds/derived/monoatomicCloud/monoatomicCloud.H
View File

@ -1,52 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::monoatomicCloud
Description
Cloud class to simulate monoatomic molecules
SourceFiles
monoatomicCloud.C
\*---------------------------------------------------------------------------*/
#ifndef monoatomicCloud_H
#define monoatomicCloud_H
#include "MoleculeCloud.H"
#include "monoatomic.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef MoleculeCloud<monoatomic> monoatomicCloud;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

52
src/lagrangian/molecularDynamics/molecule/clouds/derived/polyatomicCloud/polyatomicCloud.H
View File

@ -1,52 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::polyatomicCloud
Description
Cloud class to simulate polyatomic molecules
SourceFiles
polyatomicCloud.C
\*---------------------------------------------------------------------------*/
#ifndef polyatomicCloud_H
#define polyatomicCloud_H
#include "MoleculeCloud.H"
#include "polyatomic.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
typedef MoleculeCloud<polyatomic> polyatomicCloud;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

573
src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.C
View File

@ -1,573 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "controllers.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::controllers::controllers
(
const polyMesh& mesh
)
:
time_(mesh.time()),
controllersDict_
(
IOobject
(
"controllersDict",
time_.system(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
)
),
stateControllersList_(),
sCNames_(),
sCIds_(),
sCFixedPathNames_(),
stateControllers_(),
fluxControllersList_(),
fCNames_(),
fCIds_(),
fCFixedPathNames_(),
fluxControllers_()
{}
Foam::controllers::controllers
(
const polyMesh& mesh,
polyatomicCloud& cloud
)
:
time_(mesh.time()),
controllersDict_
(
IOobject
(
"controllersDict",
time_.system(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
stateControllersList_(controllersDict_.lookup("stateControllers")),
sCNames_(stateControllersList_.size()),
sCIds_(stateControllersList_.size()),
sCFixedPathNames_(stateControllersList_.size()),
stateControllers_(stateControllersList_.size()),
fluxControllersList_(controllersDict_.lookup("fluxControllers")),
fCNames_(fluxControllersList_.size()),
fCIds_(fluxControllersList_.size()),
fCFixedPathNames_(fluxControllersList_.size()),
fluxControllers_(fluxControllersList_.size())
{
Info << nl << "Creating controllers" << nl << endl;
// state controllers
if (!stateControllers_.empty())
{
forAll(stateControllers_, sC)
{
const entry& controllersI = stateControllersList_[sC];
const dictionary& controllersIDict = controllersI.dict();
stateControllers_[sC] = autoPtr<stateController>
(
stateController::New(time_, cloud, controllersIDict)
);
sCNames_[sC] = stateControllers_[sC]->type();
sCIds_[sC] = sC;
}
}
//- flux controllers
if (!fluxControllers_.empty())
{
forAll(fluxControllers_, fC)
{
const entry& controllersI = fluxControllersList_[fC];
const dictionary& controllersIDict = controllersI.dict();
fluxControllers_[fC] = autoPtr<fluxController>
(
fluxController::New(time_, cloud, controllersIDict)
);
fCNames_[fC] = fluxControllers_[fC]->type();
fCIds_[fC] = fC;
}
}
// creating directories for state controllers
if (!nStateControllers_.empty())
{
// case/controllers
fileName controllersPath(time_.path()/"controllers");
if (!isDir(controllersPath))
{
mkDir(controllersPath);
}
// case/controllers/<cloudName>
fileName controllersPath(controllersPath/cloud.name());
if (!isDir(controllersPath))
{
mkDir(controllersPath);
}
// case/controllers/<cloudName>/stateControllers
fileName stateControllersPath(controllersPath/"stateControllers");
if (!isDir(stateControllersPath))
{
mkDir(stateControllersPath);
}
forAll(stateControllers_, sC)
{
if (stateControllers_[sC]->writeInCase())
{
// case/controllers/<cloudName>/
// stateControllers/<stateControllerModel>
fileName stateControllerPath(stateControllersPath/sCNames_[sC]);
if (!isDir(stateControllerPath))
{
mkDir(stateControllerPath);
}
const word& regionName = stateControllers_[sC]->regionName();
// case/controllers/<cloudName>/
// stateControllers/<stateControllerModel>/<cellZoneName>
fileName zonePath(stateControllerPath/regionName);
if (!isDir(zonePath))
{
mkDir(zonePath);
}
sCFixedPathNames_[sC] = zonePath;
}
}
}
// creating directories for flux controllers
if (nFluxControllers_ > 0)
{
// case/controllers
fileName controllersPath(time_.path()/"controllers");
if ( !isDir(controllersPath) )
{
mkDir(controllersPath);
}
// case/controllers/<cloudName>
fileName controllersPath(time_.path()/cloud.name());
if ( !isDir(controllersPath) )
{
mkDir(controllersPath);
}
// case/controllers/<cloudName>/fluxControllers
fileName fluxControllersPath(controllersPath/"fluxControllers");
if (!isDir(fluxControllersPath))
{
mkDir(fluxControllersPath);
}
forAll(fluxControllers_, fC)
{
if (fluxControllers_[fC]->writeInCase())
{
// case/controllers/<cloudName>/
// fluxControllers/<fluxControllerModel>
fileName fluxControllerPath(fluxControllersPath/fCNames_[fC]);
if (!isDir(fluxControllerPath))
{
mkDir(fluxControllerPath);
}
const word& regionName = fluxControllers_[fC]->regionName();
// case/controllers/<cloudName>/
// fluxControllers/<fluxControllerModel>/<faceZoneName>
fileName zonePath(fluxControllerPath/regionName);
if (!isDir(zonePath))
{
mkDir(zonePath);
}
fCFixedPathNames_[fC] = zonePath;
}
}
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
controllers::~controllers()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void controllers::initialConfig()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->initialConfiguration();
}
forAll(fluxControllers_, fC)
{
fluxControllers_[fC]->initialConfiguration();
}
}
void controllers::updateTimeInfo()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->updateTime();
}
forAll(fluxControllers_, fC)
{
fluxControllers_[fC]->updateTime();
}
}
void controllers::controlState()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->controlMols();
}
}
void controllers::controlVelocitiesI()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->controlMolsBeg();
}
}
void controllers::controlVelocitiesII()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->controlMolsEnd();
}
}
void controllers::controlPriorToForces()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->controlBeforeForces();
}
}
void controllers::calculateStateProps()
{
forAll(stateControllers_, sC)
{
stateControllers_[sC]->calculateProperties();
}
forAll(fluxControllers_, fC)
{
fluxControllers_[fC]->calculateProperties();
}
}
void controllers::outputStateResults()
{
const Time& runTime = time_;
if (runTime.outputTime())
{
// creating a set of directories in the current time directory
{
List<fileName> timePathNames(sCFixedPathNames_.size());
if (nStateControllers_ > 0)
{
if (Pstream::master())
{
// case/<timeDir>/uniform
fileName uniformTimePath
(
runTime.path()/runTime.timeName()/"uniform"
);
if (!isDir(uniformTimePath))
{
mkDir(uniformTimePath);
}
if (!stateControllers_.empty())
{
// case/<timeDir>/uniform/controllers
fileName controllersTimePath
(
uniformTimePath/"controllers"
);
if (!isDir(controllersTimePath))
{
mkDir(controllersTimePath);
}
// case/<timeDir>/uniform/controllers/<cloudName>
fileName cloudTimePath
(
controllersTimePath/cloud.name()
);
if (!isDir(cloudTimePath))
{
mkDir(cloudTimePath);
}
// case/<timeDir>/uniform/controllers/<cloudName>/
fileName stateControllersTimePath
(
cloudTimePath/"stateControllers"
);
if (!isDir(stateControllersTimePath))
{
mkDir(stateControllersTimePath);
}
forAll(stateControllers_, sC)
{
if (stateControllers_[sC]->writeInTimeDir())
{
// case/<timeDir>/uniform/controllers/
// <cloudName>/<stateControllerModel>
fileName sCTimePath
(
stateControllersTimePath/sCNames_[sC]
);
if (!isDir(sCTimePath))
{
mkDir(sCTimePath);
}
// Creating directory for different zones but
// of the same model
const word& regionName =
stateControllers_[sC]->regionName();
// case/<timeDir>/uniform/controllers/
// <cloudName>/<stateControllerModel>/
// <cellZoneName>
fileName zoneTimePath(sCTimePath/regionName);
if (!isDir(zoneTimePath))
{
mkDir(zoneTimePath);
}
timePathNames[sC] = zoneTimePath;
}
}
}
}
}
// write out data
forAll(stateControllers_, sC)
{
stateControllers_[sC]->output
(
sCFixedPathNames_[sC],
timePathNames[sC]
);
}
}
{
List<fileName> timePathNames(fCFixedPathNames_.size());
if (nFluxControllers_ > 0)
{
if (Pstream::master())
{
// case/<timeDir>/uniform
fileName uniformTimePath
(
runTime.path()/runTime.timeName()/"uniform"
);
if (!isDir(uniformTimePath))
{
mkDir(uniformTimePath);
}
if (!fluxControllers_.empty())
{
// case/<timeDir>/uniform/controllers
fileName controllersTimePath
(
uniformTimePath/"controllers"
);
if (!isDir(controllersTimePath))
{
mkDir(controllersTimePath);
}
// case/<timeDir>/uniform/controllers/<cloudName>
fileName cloudTimePath
(
controllersTimePath/cloud.name()
);
if (!isDir(cloudTimePath))
{
mkDir(cloudTimePath);
}
// case/<timeDir>/uniform/fluxControllers
fileName controllersTimePath
(
cloudTimePath/"fluxControllers"
);
if (!isDir(controllersTimePath))
{
mkDir(controllersTimePath);
}
forAll(fluxControllers_, fC)
{
if (stateControllers_[fC]->writeInTimeDir())
{
// case/<timeDir>/uniform/controllers/
// <cloudName>/<fluxControllerModel>
fileName fCTimePath
(
controllersTimePath/fCNames_[fC]
);
if (!isDir(fCTimePath))
{
mkDir(fCTimePath);
}
const word& regionName =
fluxControllers_[fC]->regionName();
// case/<timeDir>/uniform/controllers/
// <cloudName>/<fluxControllerModel>/
// <faceZoneName>
fileName zoneTimePath(fCTimePath/regionName);
if (!isDir(zoneTimePath))
{
mkDir(zoneTimePath);
}
timePathNames[fC] = zoneTimePath;
}
}
}
}
}
// write out data
forAll(fluxControllers_, fC)
{
fluxControllers_[fC]->output
(
fCFixedPathNames_[fC],
timePathNames[fC]
);
}
}
// Re-read dictionaries for modified properties (run-time selection)
{
stateControllersList_.clear();
stateControllersList_ = controllersDict_.lookup("stateControllers");
forAll(stateControllers_, sC)
{
const entry& controllersI = stateControllersList_[sC];
const dictionary& controllersIDict = controllersI.dict();
stateControllers_[sC]->updateProperties(controllersIDict);
}
}
{
fluxControllersList_.clear();
fluxControllersList_ = controllersDict_.lookup("fluxControllers");
forAll(fluxControllers_, fC)
{
const entry& controllersI = fluxControllersList_[fC];
const dictionary& controllersIDict = controllersI.dict();
fluxControllers_[fC]->updateProperties(controllersIDict);
}
}
}
}
// ************************************************************************* //

163
src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllers.H
View File

@ -1,163 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
controllers
Description
Stores all the information for the controllers models defined within
the controllersDict, and selects & builds the models automatically.
\*---------------------------------------------------------------------------*/
#ifndef controllers_H
#define controllers_H
#include "List.H"
#include "IOdictionary.H"
#include "autoPtr.H"
#include "polyMesh.H"
#include "timeData.H"
#include "stateController.H"
#include "fluxController.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class controllers Declaration
\*---------------------------------------------------------------------------*/
class controllers
{
// Private data
Time& time_;
//- The entire dictionary (containing multiple subDictionaries)
IOdictionary controllersDict_;
//- state controllers
PtrList<entry> stateControllersList_;
List<word> sCNames_;
List<label> sCIds_;
List<fileName> sCFixedPathNames_;
List< autoPtr<stateController> > stateControllers_;
//- flux controllers
PtrList<entry> fluxControllersList_;
List<word> fCNames_;
List<label> fCIds_;
List<fileName> fCFixedPathNames_;
List< autoPtr<fluxController> > fluxControllers_;
public:
// Constructors
//- Null Constructor
controllers
(
const polyMesh& mesh
);
//- Constructor for with cloud
controllers
(
const polyMesh& mesh,
polyatomicCloud& cloud
);
//- Destructor
~controllers();
// Member Functions
//- Initial configuration call this function after the polyatomicCloud
// is completely initialised
void initialConfig();
//- this function is to be called at the beginning of the MD time-step.
// since we have placed a non-referenced time-data class in the
// state-controller class.
void updateTimeInfo();
//- control molecular state -- call this after the intermolecular force
// calulation
void controlState();
//-
void controlVelocitiesI();
//-
void controlVelocitiesII();
//-
void controlPriorToForces();
//- calculate properties -- call this at the end of the MD time-step.
void calculateStateProps();
//- output -- call this function at the end of the MD time-step
void outputStateResults();
// Access
//-
inline List< autoPtr<stateController> >& stateControllers();
//-
inline const List< autoPtr<stateController> >&
stateControllers() const;
//-
inline List< autoPtr<fluxController> >& fluxControllers();
//-
inline const List< autoPtr<fluxController> >&
fluxControllers() const;
//-
inline const List<word>& stateControllersNames() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "controllersI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

64
src/lagrangian/molecularDynamics/molecule/controllers/basic/controllers/controllersI.H
View File

@ -1,64 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::List<Foam::autoPtr<Foam::stateController> >&
Foam::controllers::stateControllers()
{
return stateControllers_;
}
const Foam::List<Foam::autoPtr<Foam::stateController> >&
Foam::controllers::stateControllers() const
{
return stateControllers_;
}
Foam::List<Foam::autoPtr<Foam::fluxController> >&
Foam::controllers::fluxControllers()
{
return fluxControllers_;
}
const Foam::List< autoPtr<fluxController> >&
Foam::controllers::fluxControllers() const
{
return fluxControllers_;
}
const Foam::List<Foam::word>& Foam::controllers::stateControllersNames() const
{
return sCNames_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

524
src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.C
View File

@ -1,524 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "waterFluxController.H"
#include "IFstream.H"
#include "graph.H"
#include "polyatomicCloud.H"
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(waterFluxController, 0);
defineRunTimeSelectionTable(waterFluxController, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
waterFluxController::waterFluxController
(
Time& t,
polyatomicCloud& cloud,
const dictionary& dict
)
:
mesh_(refCast<const fvMesh>(cloud.mesh())),
cloud_(cloud),
rndGen_(clock::getTime()),
controllerDict_(dict.subDict("controllerProperties")),
timeDict_(controllerDict_.subDict("timeProperties")),
time_(t, timeDict_),
regionName_(controllerDict_.lookup("zoneName")),
regionId_(-1),
zoneSurfaceArea_(0.0),
internalFaces_(),
processorFaces_(),
control_(true),
readStateFromFile_(true),
singleValueController_(false),
density_(0.0),
velocity_(vector::zero),
temperature_(0.0),
pressure_(0.0),
strainRate_(tensor::zero),
tempGradient_(vector::zero),
fieldController_(false),
densities_(),
velocities_(),
temperatures_(),
pressures_(),
writeInTimeDir_(true),
writeInCase_(true)
{
const faceZoneMesh& faceZones = mesh_.faceZones();
regionId_ = faceZones.findZoneID(regionName_);
if (regionId_ == -1)
{
FatalErrorIn("waterFluxController::waterFluxController()")
<< "Cannot find region (faceZone): " << regionName_ << nl << "in: "
<< t.time().system()/"controllersDict"
<< exit(FatalError);
}
control_ = Switch(controllerDict_.lookup("controlSwitch"));
readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
setFacesInfo();
}
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
autoPtr<waterFluxController> waterFluxController::New
(
Time& t,
polyatomicCloud& cloud,
const dictionary& dict
)
{
word waterFluxControllerName
(
dict.lookup("fluxControllerModel")
);
Info<< "Selecting fluxController "
<< waterFluxControllerName << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(waterFluxControllerName);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalError
<< "waterFluxController::New(const dictionary&) : " << endl
<< " unknown waterFluxController type "
<< waterFluxControllerName
<< ", constructor not in hash table" << endl << endl
<< " Valid injector types are :" << endl;
Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
}
return autoPtr<waterFluxController>
(
cstrIter()(t, cloud, dict)
);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
waterFluxController::~waterFluxController()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// void waterFluxController::updateTime()
// {
// time_++;
//
// const scalar& t = time_.time().timeOutputValue();
//
// if ((t - initialTime_) < timePeriod_)
// {
// time_.controlTimeInterval().endTime() = false;
// // control_ = false;
// }
// else
// {
// // control_ = true;
// }
// }
void waterFluxController::setFacesInfo()
{
const labelList& faces = controlZone();
if (Pstream::parRun())
{
DynamicList<label> processorFaces(0);
forAll(mesh_.boundaryMesh(), patchI)
{
const polyPatch& patch = mesh_.boundaryMesh()[patchI];
if (isA<processorPolyPatch>(patch))
{
for (label p = 0; p < patch.size(); p++)
{
label patchFaceI = p + patch.start();
label faceId = findIndex (faces, patchFaceI);
if (faceId != -1)
{
processorFaces.append(patchFaceI);
}
}
}
}
processorFaces.shrink();
processorFaces_.setSize(processorFaces.size(), -1);
forAll(processorFaces, f)
{
processorFaces_[f] = processorFaces[f];
}
label nInternalFaces = faces.size() - processorFaces.size();
internalFaces_.setSize(nInternalFaces, -1);
label counter = 0;
forAll(faces, f)
{
const label& faceI = faces[f];
if (findIndex(processorFaces, faceI) == -1)
{
internalFaces_[counter] = faceI;
counter++;
}
}
// Pout << "processorFaces: " << processorFaces_ << endl;
// Pout << "internalFaces: " << internalFaces_ << endl;
forAll(internalFaces_, f)
{
const label& faceI = internalFaces_[f];
zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
}
// faces on a zone located on a processor cut belong to both processors
// (hence the 0.5)
forAll(processorFaces_, f)
{
const label& faceI = processorFaces_[f];
zoneSurfaceArea_ += 0.5*mag(mesh_.faceAreas()[faceI]);
}
if (Pstream::parRun())
{
for (int p = 0; p < Pstream::nProcs(); p++)
{
if (p != Pstream::myProcNo())
{
const int proc = p;
{
OPstream toNeighbour(Pstream::blocking, proc);
toNeighbour << zoneSurfaceArea_;
}
}
}
//- receiving
for (int p = 0; p < Pstream::nProcs(); p++)
{
if (p != Pstream::myProcNo())
{
scalar zoneSurfaceAreaProc;
const int proc = p;
{
IPstream fromNeighbour(Pstream::blocking, proc);
fromNeighbour >> zoneSurfaceAreaProc;
}
zoneSurfaceArea_ += zoneSurfaceAreaProc;
}
}
}
}
else
{
forAll(faces, f)
{
const label& faceI = faces[f];
zoneSurfaceArea_ += mag(mesh_.faceAreas()[faceI]);
}
}
}
void waterFluxController::updateTime()
{
time_++;
// const scalar& t = time_.time().timeOutputValue();
//
// if ((t - initialTime_) < timePeriod_)
// {
// time_.controlTimeInterval().endTime() = false;
// // control_ = false;
// }
// else
// {
// // control_ = true;
// }
}
void waterFluxController::updateFluxControllerProperties
(
const dictionary& newDict
)
{
controllerDict_ = newDict.subDict("controllerProperties");
//- you can reset the controlling zone from here. This essentially
// means that the coupling zone can infact move arbitrarily. To make
// this happen we probably need to devise a technique for automatically
// changing the cellZone else where, and then calling this function to
// reset the controlling zone in which the controller operates in.
if (controllerDict_.found("controlSwitch"))
{
control_ = Switch(controllerDict_.lookup("controlSwitch"));
}
if (controllerDict_.found("readStateFromFile"))
{
readStateFromFile_ = Switch
(
controllerDict_.lookup("readStateFromFile")
);
}
}
const labelList& waterFluxController::controlZone() const
{
return mesh_.faceZones()[regionId_];
}
label waterFluxController::isFaceOnControlZone(const label& faceI)
{
const label f = findIndex(controlZone(), faceI);
return f;
}
const word& waterFluxController::regionName() const
{
return regionName_;
}
const scalar& waterFluxController::density() const
{
return density_;
}
scalar& waterFluxController::density()
{
return density_;
}
const vector& waterFluxController::velocity() const
{
return velocity_;
}
vector& waterFluxController::velocity()
{
return velocity_;
}
const scalar& waterFluxController::temperature() const
{
return temperature_;
}
scalar& waterFluxController::temperature()
{
return temperature_;
}
const scalar& waterFluxController::pressure() const
{
return pressure_;
}
scalar& waterFluxController::pressure()
{
return pressure_;
}
const tensor& waterFluxController::strainRate() const
{
return strainRate_;
}
tensor& waterFluxController::strainRate()
{
return strainRate_;
}
const vector& waterFluxController::tempGradient() const
{
return tempGradient_;
}
vector& waterFluxController::tempGradient()
{
return tempGradient_;
}
const scalarField& waterFluxController::densityField() const
{
return densities_;
}
scalarField& waterFluxController::densityField()
{
return densities_;
}
const vectorField& waterFluxController::velocityField() const
{
return velocities_;
}
vectorField& waterFluxController::velocityField()
{
return velocities_;
}
const scalarField& waterFluxController::temperatureField() const
{
return temperatures_;
}
scalarField& waterFluxController::temperatureField()
{
return temperatures_;
}
const scalarField& waterFluxController::pressureField() const
{
return pressures_;
}
scalarField& waterFluxController::pressureField()
{
return pressures_;
}
const bool& waterFluxController::singleValueController() const
{
return singleValueController_;
}
bool& waterFluxController::singleValueController()
{
return singleValueController_;
}
const bool& waterFluxController::fieldController() const
{
return fieldController_;
}
bool& waterFluxController::fieldController()
{
return fieldController_;
}
const bool& waterFluxController::writeInTimeDir() const
{
return writeInTimeDir_;
}
const bool& waterFluxController::writeInCase() const
{
return writeInCase_;
}
// const scalar waterFluxController::avReqDensity() const
// {
// scalar totalDensity = 0.0;
//
// forAll(densities_, c)
// {
// totalDensity += densities_[c];
// }
//
// if (cells_.size() > 0)
// {
// totalDensity /= scalar(cells_.size());
// }
//
// return totalDensity;
// }
//
// const vector waterFluxController::avReqVelocity() const
// {
// vector totalVel = vector::zero;
//
// forAll(velocities_, c)
// {
// totalVel += velocities_[c];
// }
//
// if (cells_.size() > 0)
// {
// totalVel /= scalar(cells_.size());
// }
//
// return totalVel;
// }
//
// const scalar waterFluxController::avReqTemperature() const
// {
// scalar totalTemp = 0.0;
//
// forAll(densities_, c)
// {
// totalTemp += temperatures_[c];
// }
//
// if (cells_.size() > 0)
// {
// totalTemp /= scalar(cells_.size());
// }
//
// return totalTemp;
// }
} // End namespace Foam
// ************************************************************************* //

272
src/lagrangian/molecularDynamics/molecule/controllers/basic/fluxController/fluxController.H
View File

@ -1,272 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
waterFluxController
Description
SourceFiles
waterFluxControllerI.H
waterFluxController.C
waterFluxControllerIO.C
\*---------------------------------------------------------------------------*/
#ifndef waterFluxController_H
#define waterFluxController_H
#include "IOdictionary.H"
#include "Time.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "vector.H"
#include "volFields.H"
#include "Random.H"
#include "polyatomic.H"
#include "timeData.H"
#include "writeTimeData.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class waterFluxController Declaration
\*---------------------------------------------------------------------------*/
class waterFluxController
{
protected:
// Protected data
// Time& time_;
const fvMesh& mesh_;
polyatomicCloud& cloud_;
Random rndGen_;
//- subDictionary containing the properties
dictionary controllerDict_;
dictionary timeDict_;
timeData time_;
//- name of face zone
word regionName_;
label regionId_;
// labelList faces_;
scalar zoneSurfaceArea_;
labelList internalFaces_;
labelList processorFaces_;
bool control_;
bool readStateFromFile_;
//- set all the properties below from model if required
bool singleValueController_;
// target values
scalar density_;
vector velocity_;
scalar temperature_;
scalar pressure_;
tensor strainRate_;
vector tempGradient_;
bool fieldController_;
//- targeted fields
scalarField densities_;
vectorField velocities_;
scalarField temperatures_;
scalarField pressures_;
bool writeInTimeDir_;
bool writeInCase_;
// Private Member Functions
void setFacesInfo();
public:
//- Runtime type information
TypeName("waterFluxController");
// Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
waterFluxController,
dictionary,
(
Time& t,
polyatomicCloud& cloud,
const dictionary& dict
),
(t, cloud, dict)
);
// Constructors
//- Construct from components
waterFluxController
(
Time& t,
polyatomicCloud& cloud,
const dictionary& dict
);
// Selectors
static autoPtr<waterFluxController> New
(
Time& t,
polyatomicCloud& cloud,
const dictionary& dict
);
// Destructor
virtual ~waterFluxController();
// Member Functions
void updateTime();
//- create an initial configuration
virtual void initialConfiguration() = 0;
//- calculate any required properties
virtual void calculateProperties() = 0;
//- control the polyatomic from the tracking function
virtual void controlMol
(
polyatomic& mol,
polyatomic::trackData& td
) = 0;
//- output data
virtual void output
(
const fileName& fixedPathName,
const fileName& timePath
) = 0;
//- E. update properties from a modified dictionary
virtual void updateProperties(const dictionary&) = 0;
void updateFluxControllerProperties(const dictionary&);
// Access
//- return the control zone cells
const labelList& controlZone() const;
label isFaceOnControlZone(const label& faceI);
//- return the control zone name
const word& regionName() const;
//- return the targeted values
const scalar& density() const;
scalar& density();
const vector& velocity() const;
vector& velocity();
const scalar& temperature() const;
scalar& temperature();
const scalar& pressure() const;
scalar& pressure();
const tensor& strainRate() const;
tensor& strainRate();
const vector& tempGradient() const;
vector& tempGradient();
//- return the targeted fields
const scalarField& densityField() const;
scalarField& densityField();
const vectorField& velocityField() const;
vectorField& velocityField();
const scalarField& temperatureField() const;
scalarField& temperatureField();
const scalarField& pressureField() const;
scalarField& pressureField();
const bool& singleValueController() const;
bool& singleValueController();
const bool& fieldController() const;
bool& fieldController();
const bool& writeInTimeDir() const;
const bool& writeInCase() const;
// const scalar avReqDensity() const;
// const vector avReqVelocity() const;
// const scalar avReqTemperature() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

490
src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.C
View File

@ -1,490 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "stateController.H"
#include "IFstream.H"
#include "polyatomicCloud.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(stateController, 0);
defineRunTimeSelectionTable(stateController, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::stateController::stateController
(
polyatomicCloud& cloud,
const dictionary& dict
)
:
mesh_(refCast<const fvMesh>(cloud.mesh())),
cloud_(cloud),
rndGen_(clock::getTime()),
controllerDict_(dict.subDict("controllerProperties")),
timeDict_(controllerDict_.subDict("timeProperties")),
time_(mesh_.time(), timeDict_),
timePeriod_(readScalar(timeDict_.lookup("initialTimePeriod"))), //temp
initialTime_(time_.time().startTime().value()),
regionName_(controllerDict_.lookup("zoneName")),
regionId_(-1),
control_(true),
readStateFromFile_(true),
singleValueController_(false),
density_(0.0),
velocity_(vector::zero),
temperature_(0.0),
pressure_(0.0),
strainRate_(tensor::zero),
tempGradient_(vector::zero),
fieldController_(false),
densities_(),
velocities_(),
temperatures_(),
pressures_(),
writeInTimeDir_(true),
writeInCase_(true)
{
const cellZoneMesh& cellZones = mesh_.cellZones();
regionId_ = cellZones.findZoneID(regionName_);
if (regionId_ == -1)
{
FatalErrorIn("stateController::stateController()")
<< "Cannot find region: " << regionName_ << nl << "in: "
<< time_.time().system()/"controllersDict"
<< exit(FatalError);
}
control_ = Switch(controllerDict_.lookup("controlSwitch"));
readStateFromFile_ = Switch(controllerDict_.lookup("readStateFromFile"));
const scalar& avTimeInterval = time_.averageTimeInterval().deltaT();
if ((timePeriod_ < avTimeInterval) && (timePeriod_ > 0.0))
{
timePeriod_ = avTimeInterval;
}
}
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::stateController> Foam::stateController::New
(
polyatomicCloud& cloud,
const dictionary& dict
)
{
word stateControllerName
(
dict.lookup("stateControllerModel")
);
Info<< "Selecting stateController "
<< stateControllerName << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(stateControllerName);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalError
<< "stateController::New(const dictionary&) : " << endl
<< " unknown stateController type "
<< stateControllerName
<< ", constructor not in hash table" << endl << endl
<< " Valid types are :" << endl;
Info<< dictionaryConstructorTablePtr_->toc() << abort(FatalError);
}
return autoPtr<stateController>
(
cstrIter()(cloud, dict)
);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::stateController::~stateController()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::stateController::updateTime()
{
time_++;
const scalar& t = time_.time().timeOutputValue();
if ((t - initialTime_) < timePeriod_)
{
time_.controlTimeInterval().endTime() = false;
}
}
void Foam::stateController::updateStateControllerProperties
(
const dictionary& newDict
)
{
controllerDict_ = newDict.subDict("controllerProperties");
if (controllerDict_.found("controlSwitch"))
{
control_ = Switch(controllerDict_.lookup("controlSwitch"));
}
if (controllerDict_.found("readStateFromFile"))
{
readStateFromFile_ = Switch
(
controllerDict_.lookup("readStateFromFile")
);
}
timeDict_ = controllerDict_.subDict("timeProperties");
if (timeDict_.found("resetAtOutput"))
{
time_.resetFieldsAtOutput() = Switch(timeDict_.lookup("resetAtOutput"));
}
}
const Foam::labelList& Foam::stateController::controlZone() const
{
return mesh_.cellZones()[regionId_];
}
const Foam::word& Foam::stateController::regionName() const
{
return regionName_;
}
Foam::scalar Foam::stateController::density() const
{
return density_;
}
Foam::scalar& Foam::stateController::density()
{
return density_;
}
const Foam::vector& Foam::stateController::velocity() const
{
return velocity_;
}
Foam::vector& Foam::stateController::velocity()
{
return velocity_;
}
Foam::scalar Foam::stateController::temperature() const
{
return temperature_;
}
Foam::scalar& Foam::stateController::temperature()
{
return temperature_;
}
const Foam::scalar& Foam::stateController::pressure() const
{
return pressure_;
}
Foam::scalar& Foam::stateController::pressure()
{
return pressure_;
}
const Foam::tensor& Foam::stateController::strainRate() const
{
return strainRate_;
}
Foam::tensor& Foam::stateController::strainRate()
{
return strainRate_;
}
const Foam::vector& Foam::stateController::tempGradient() const
{
return tempGradient_;
}
Foam::vector& Foam::stateController::tempGradient()
{
return tempGradient_;
}
const Foam::scalarField& Foam::stateController::densityField() const
{
return densities_;
}
Foam::scalarField& Foam::stateController::densityField()
{
return densities_;
}
const Foam::vectorField& Foam::stateController::velocityField() const
{
return velocities_;
}
Foam::vectorField& Foam::stateController::velocityField()
{
return velocities_;
}
const Foam::scalarField& Foam::stateController::temperatureField() const
{
return temperatures_;
}
Foam::scalarField& Foam::stateController::temperatureField()
{
return temperatures_;
}
const Foam::scalarField& Foam::stateController::pressureField() const
{
return pressures_;
}
Foam::scalarField& Foam::stateController::pressureField()
{
return pressures_;
}
bool Foam::stateController::singleValueController() const
{
return singleValueController_;
}
bool& Foam::stateController::singleValueController()
{
return singleValueController_;
}
bool Foam::stateController::fieldController() const
{
return fieldController_;
}
bool& Foam::stateController::fieldController()
{
return fieldController_;
}
bool Foam::stateController::writeInTimeDir() const
{
return writeInTimeDir_;
}
bool Foam::stateController::writeInCase() const
{
return writeInCase_;
}
Foam::scalar Foam::stateController::avReqDensity()
{
scalar totalDensity = 0.0;
if (singleValueController_)
{
totalDensity = density_;
}
else if (fieldController_)
{
label controlCells = controlZone().size();
forAll(densities_, c)
{
totalDensity += densities_[c];
}
if (Pstream::parRun())
{
reduce(totalDensity, sumOp<scalar>());
reduce(controlCells, sumOp<label>());
}
if (controlCells > 0)
{
totalDensity /= scalar(controlCells);
}
}
return totalDensity;
}
Foam::vector Foam::stateController::avReqVelocity()
{
vector totalVel = vector::zero;
if (singleValueController_)
{
totalVel = velocity_;
}
else if (fieldController_)
{
label controlCells = controlZone().size();
forAll(velocities_, c)
{
totalVel += velocities_[c];
}
if (Pstream::parRun())
{
reduce(totalVel, sumOp<vector>());
reduce(controlCells, sumOp<label>());
}
if (controlCells > 0)
{
totalVel /= scalar(controlCells);
}
}
return totalVel;
}
Foam::scalar Foam::stateController::avReqTemperature()
{
scalar totalTemp = 0.0;
if (singleValueController_)
{
totalTemp = temperature_;
}
else if (fieldController_)
{
label controlCells = controlZone().size();
forAll(temperatures_, c)
{
totalTemp += temperatures_[c];
}
if (Pstream::parRun())
{
reduce(totalTemp, sumOp<scalar>());
reduce(controlCells, sumOp<label>());
}
if (controlCells > 0)
{
totalTemp /= scalar(controlCells);
}
}
return totalTemp;
}
Foam::scalar Foam::stateController::avReqPressure()
{
scalar totalPressure = 0.0;
if (singleValueController_)
{
totalPressure = pressure_;
}
else if (fieldController_)
{
label controlCells = controlZone().size();
forAll(pressures_, c)
{
totalPressure += pressures_[c];
}
if (Pstream::parRun())
{
reduce(totalPressure, sumOp<scalar>());
reduce(controlCells, sumOp<label>());
}
if (controlCells > 0)
{
totalPressure /= scalar(controlCells);
}
}
return totalPressure;
}
// ************************************************************************* //

270
src/lagrangian/molecularDynamics/molecule/controllers/basic/stateController/stateController.H
View File

@ -1,270 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
stateController
Description
Basic/abstract class of a state controller
SourceFiles
stateControllerI.H
stateController.C
stateControllerIO.C
\*---------------------------------------------------------------------------*/
#ifndef stateController_H
#define stateController_H
#include "IOdictionary.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
#include "vector.H"
#include "volFields.H"
#include "Random.H"
#include "polyatomic.H"
#include "timeData.H"
#include "writeTimeData.H"
#include "selectIds.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class stateController Declaration
\*---------------------------------------------------------------------------*/
class stateController
{
protected:
// Protected data
//-
const fvMesh& mesh_;
//-
polyatomicCloud& cloud_;
//-
Random rndGen_;
//- subDictionary containing the properties
dictionary controllerDict_;
//-
dictionary timeDict_;
//-
timeData time_;
//-
scalar timePeriod_;
//-
scalar initialTime_;
//- name of control zone
word regionName_;
//-
label regionId_;
//-
bool control_;
//-
bool readStateFromFile_;
//- set all the properties below from model if required
//-
bool singleValueController_;
//- target values
scalar density_;
vector velocity_;
scalar temperature_;
scalar pressure_;
tensor strainRate_;
vector tempGradient_;
//- set this in model
bool fieldController_;
//- targeted fields
scalarField densities_;
vectorField velocities_;
scalarField temperatures_;
scalarField pressures_;
// set these in model
bool writeInTimeDir_;
bool writeInCase_;
public:
//- Runtime type information
TypeName("stateController");
//- Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
stateController,
dictionary,
(
polyatomicCloud& cloud,
const dictionary& dict
),
(t, cloud, dict)
);
// Constructors
//- Construct from components
stateController
(
polyatomicCloud& cloud,
const dictionary& dict
);
// Selectors
static autoPtr<stateController> New
(
polyatomicCloud& cloud,
const dictionary& dict
);
// Destructor
virtual ~stateController();
// Member Functions
void updateTime();
//- create an initial configuration
virtual void initialConfiguration() = 0;
//- calculate any required properties
virtual void calculateProperties() = 0;
//- control molecules at different stages of the integration time-step
virtual void controlMolsBeg() = 0;
virtual void controlBeforeForces() = 0;
virtual void controlMols() = 0;
virtual void controlMolsEnd() = 0;
//- output data
virtual void output
(
const fileName& fixedPathName,
const fileName& timePath
) = 0;
//- update properties from a modified dictionary
virtual void updateProperties(const dictionary&) = 0;
void updateStateControllerProperties(const dictionary&);
// Access
//- return the control zone cells
const labelList& controlZone() const;
//- return the control zone name
const word& regionName() const;
//- return the targeted fields
scalar density() const;
scalar& density();
const vector& velocity() const;
vector& velocity();
scalar temperature() const;
scalar& temperature();
scalar pressure() const;
scalar& pressure();
const tensor& strainRate() const;
tensor& strainRate();
const vector& tempGradient() const;
vector& tempGradient();
//- return the targeted fields
const scalarField& densityField() const;
scalarField& densityField();
const vectorField& velocityField() const;
vectorField& velocityField();
const scalarField& temperatureField() const;
scalarField& temperatureField();
const scalarField& pressureField() const;
scalarField& pressureField();
bool singleValueController() const;
bool& singleValueController();
bool fieldController() const;
bool& fieldController();
bool writeInTimeDir() const;
bool writeInCase() const;
scalar avReqDensity();
vector avReqVelocity();
scalar avReqTemperature();
scalar avReqPressure();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

114
src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
View File

@ -1,114 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constPropSite.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::constPropSite::constPropSite()
:
siteReferencePosition_(vector::zero),
siteMass_(0.0),
siteCharge_(0.0),
siteId_(0),
name_(),
pairPotentialSite_(false),
electrostaticSite_(false)
{}
Foam::constPropSite::constPropSite
(
const vector& siteReferencePosition,
const scalar& siteMass,
const scalar& siteCharge,
const label& siteId,
const word& name,
const bool& pairPotentialSite,
const bool& electrostaticSite
)
:
siteReferencePosition_(siteReferencePosition),
siteMass_(siteMass),
siteCharge_(siteCharge),
siteId_(siteId),
name_(name),
pairPotentialSite_(pairPotentialSite),
electrostaticSite_(electrostaticSite)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::constPropSite::~constPropSite()
{}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>(Istream& is, constPropSite& cPS)
{
is >> cPS.siteReferencePosition_
>> cPS.siteMass_
>> cPS.siteCharge_
>> cPS.siteId_
>> cPS.name_
>> cPS.pairPotentialSite_
>> cPS.electrostaticSite_;
// Check state of Istream
is.check
(
"Foam::Istream& Foam::operator>>"
"(Foam::Istream&, Foam::constPropSite&)"
);
return is;
}
Foam::Ostream& Foam::operator<<(Ostream& os, const constPropSite& cPS)
{
os << token::SPACE << cPS.siteReferencePosition()
<< token::SPACE << cPS.siteMass()
<< token::SPACE << cPS.siteCharge()
<< token::SPACE << cPS.siteId()
<< token::SPACE << cPS.name()
<< token::SPACE << cPS.pairPotentialSite()
<< token::SPACE << cPS.electrostaticSite();
// Check state of Ostream
os.check
(
"Foam::Ostream& Foam::operator<<(Foam::Ostream&, "
"const Foam::constPropSite&)"
);
return os;
}
// ************************************************************************* //

185
src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
View File

@ -1,185 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Description
SourceFiles
constPropSiteI.H
constPropSite.C
\*---------------------------------------------------------------------------*/
#ifndef constPropSite_H
#define constPropSite_H
#include "vector.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class Istream;
class Ostream;
// Forward declaration of friend functions and operators
class constPropSite;
Istream& operator>>(Istream&, constPropSite&);
Ostream& operator<<(Ostream&, const constPropSite&);
/*---------------------------------------------------------------------------*\
Class constPropSite Declaration
\*---------------------------------------------------------------------------*/
class constPropSite
{
// Private data
//-
vector siteReferencePosition_;
//-
scalar siteMass_;
//-
scalar siteCharge_;
//-
label siteId_;
//-
word name_;
//-
bool pairPotentialSite_;
//-
bool electrostaticSite_;
public:
// Constructors
//- Construct null
constPropSite();
//- Construct from components
constPropSite
(
const vector& siteReferencePosition,
const scalar& siteMass,
const scalar& siteCharge,
const label& siteId,
const word& name,
const bool& pairPotentialSite,
const bool& electrostaticSite
);
// Selectors
// Destructor
~constPropSite();
// Member Functions
// Access
//-
inline const vector& siteReferencePosition() const;
//-
inline vector& siteReferencePosition();
//-
inline const scalar& siteMass() const;
//-
inline scalar& siteMass();
//-
inline const scalar& siteCharge() const;
//-
inline scalar& siteCharge();
//-
inline const label& siteId() const;
//-
inline label& siteId();
//-
inline const word& name() const;
//-
inline word& name();
//-
inline const bool& pairPotentialSite() const;
//-
inline bool& pairPotentialSite();
//-
inline const bool& electrostaticSite() const;
//-
inline bool& electrostaticSite();
// Member Operators
inline bool operator==(const constPropSite&) const;
inline bool operator!=(const constPropSite&) const;
// IOstream Operators
friend Istream& operator>>(Istream&, constPropSite&);
friend Ostream& operator<<(Ostream&, const constPropSite&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "constPropSiteI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

139
src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSiteI.H
View File

@ -1,139 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::vector& Foam::constPropSite::siteReferencePosition() const
{
return siteReferencePosition_;
}
inline Foam::vector& Foam::constPropSite::siteReferencePosition()
{
return siteReferencePosition_;
}
inline const Foam::scalar& Foam::constPropSite::siteMass() const
{
return siteMass_;
}
inline Foam::scalar& Foam::constPropSite::siteMass()
{
return siteMass_;
}
inline const Foam::scalar& Foam::constPropSite::siteCharge() const
{
return siteCharge_;
}
inline Foam::scalar& Foam::constPropSite::siteCharge()
{
return siteCharge_;
}
inline const Foam::label& Foam::constPropSite::siteId() const
{
return siteId_;
}
inline Foam::label& Foam::constPropSite::siteId()
{
return siteId_;
}
inline const Foam::word& Foam::constPropSite::name() const
{
return name_;
}
inline Foam::word& Foam::constPropSite::name()
{
return name_;
}
inline const bool& Foam::constPropSite::pairPotentialSite() const
{
return pairPotentialSite_;
}
inline bool& Foam::constPropSite::pairPotentialSite()
{
return pairPotentialSite_;
}
inline const bool& Foam::constPropSite::electrostaticSite() const
{
return electrostaticSite_;
}
inline bool& Foam::constPropSite::electrostaticSite()
{
return electrostaticSite_;
}
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
bool Foam::constPropSite::operator==
(
const constPropSite& rhs
) const
{
return
siteReferencePosition_ == rhs.siteReferencePosition_
&& siteMass_ == rhs.siteMass_
&& siteCharge_ == rhs.siteCharge_
&& siteId_ == rhs.siteId_
&& name_ == rhs.name_
&& pairPotentialSite_ == rhs.pairPotentialSite_
&& electrostaticSite_ == rhs.electrostaticSite_;
}
bool Foam::constPropSite::operator!=
(
const constPropSite& rhs
) const
{
return !(*this == rhs);
}
// ************************************************************************* //

199
src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.C
View File

@ -1,199 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*----------------------------------------------------------------------------*/
#include "monoatomic.H"
#include "Random.H"
#include "Time.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(monoatomic, 0);
defineTemplateTypeNameAndDebug(Cloud<monoatomic>, 0);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::monoatomic::move
(
monoatomic::trackingData& td,
const scalar trackTime
)
{
td.switchProcessor = false;
td.keepParticle = true;
if (special_ != SPECIAL_FROZEN)
{
return td.keepParticle;
}
const constantProperties& constProps(td.cloud().constProps(id_));
if (td.part() == trackingData::tpFirstVelocityHalfStep)
{
// First leapfrog velocity adjust part, required before tracking+force
// part
v_ += 0.5*trackTime*a_;
}
else if (td.part() == trackingData::tpLinearTrack)
{
// Leapfrog tracking part
scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{
// set the lagrangian time-step
scalar dt = min(dtMax, tEnd);
dt *= trackToFace(position() + dt*v_, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
}
setSitePositions(constProps);
}
else if (td.part() == trackingData::tpSecondVelocityHalfStep)
{
// Second leapfrog velocity adjust part, required after tracking+force
// part
a_ = siteForces_[0]/constProps.mass();
v_ += 0.5*trackTime*a_;
}
else if (td.part() != trackingData::tpRotationalTrack)
{
FatalErrorIn("monoatomic::move(trackingData&, const scalar)") << nl
<< td.part() << " is an invalid part of the integration method."
<< abort(FatalError);
}
return td.keepParticle;
}
void Foam::monoatomic::transformProperties(const tensor& T)
{
particle::transformProperties(T);
v_ = transform(T, v_);
a_ = transform(T, a_);
rf_ = transform(T, rf_);
sitePositions_[0] = position_ + (T & (sitePositions_[0] - position_));
siteForces_[0] = T & siteForces_[0];
}
void Foam::monoatomic::transformProperties(const vector& separation)
{
particle::transformProperties(separation);
if (special_ == SPECIAL_TETHERED)
{
specialPosition_ += separation;
}
sitePositions_[0] += separation;
}
void Foam::monoatomic::setSitePositions(const constantProperties& constProps)
{
sitePositions_[0] = position_;
}
void Foam::monoatomic::setSiteSizes(label size)
{
// Nothing required, size controlled internally
}
bool Foam::monoatomic::hitPatch
(
const polyPatch&,
trackingData&,
const label,
const scalar,
const tetIndices&
)
{
return false;
}
void Foam::monoatomic::hitProcessorPatch
(
const processorPolyPatch&,
trackingData& td
)
{
td.switchProcessor = true;
}
void Foam::monoatomic::hitWallPatch
(
const wallPolyPatch& wpp,
trackingData& td,
const tetIndices& tetIs
)
{
// Use of the normal from tetIs is not required as
// hasWallImpactDistance is false.
vector nw = normal();
nw /= mag(nw);
scalar vn = v_ & nw;
// Specular reflection
if (vn > 0)
{
v_ -= 2*vn*nw;
}
}
void Foam::monoatomic::hitPatch
(
const polyPatch&,
trackingData& td
)
{
td.keepParticle = false;
}
// ************************************************************************* //

418
src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomic.H
View File

@ -1,418 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::monoatomic
Description
Foam::monoatomic
SourceFiles
monoatomicI.H
monoatomic.C
monoatomicIO.C
\*---------------------------------------------------------------------------*/
#ifndef monoatomic_H
#define monoatomic_H
#include "particle.H"
#include "IOstream.H"
#include "autoPtr.H"
#include "diagTensor.H"
#include "constPropSite.H"
#include "MoleculeCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class monoatomic Declaration
\*---------------------------------------------------------------------------*/
class monoatomic
:
public particle
{
public:
// Values of special that are less than zero are for built-in functionality.
// Values greater than zero are user specifiable/expandable (i.e. test
// special_ >= SPECIAL_USER)
enum specialTypes
{
SPECIAL_TETHERED = -1,
SPECIAL_FROZEN = -2,
NOT_SPECIAL = 0,
SPECIAL_USER = 1
};
//- Class to hold monoatomic constant properties
class constantProperties
{
// Private data
//- Sites of mass, charge or interaction
FixedList<constPropSite, 1> sites_;
//- Which sites require electrostatic interactions
FixedList<label, 1> electrostaticSites_;
//- Which sites require pair interactions
FixedList<label, 1> pairPotSites_;
//-
scalar mass_;
public:
//-
inline constantProperties();
//- Construct from dictionary
inline constantProperties
(
const dictionary& dict,
const List<label>& siteIds
);
// Member functions
//-
inline const FixedList<constPropSite, 1>& sites() const;
//-
inline const FixedList<label, 1>& pairPotSites() const;
//-
inline const FixedList<label, 1>& electrostaticSites() const;
//-
inline label degreesOfFreedom() const;
//-
inline scalar mass() const;
//-
inline label nSites() const;
};
//- Class used to pass tracking data to the trackToFace function
class trackingData
:
public particle::TrackingData<MoleculeCloud<monoatomic> >
{
public:
enum trackPart
{
tpFirstVelocityHalfStep,
tpLinearTrack,
tpRotationalTrack,
tpSecondVelocityHalfStep,
tpAccess
};
private:
// label specifying which part of the integration algorithm is taking
label part_;
public:
// Constructors
trackingData(MoleculeCloud<monoatomic>& cloud, trackPart part)
:
particle::TrackingData<MoleculeCloud<monoatomic> >(cloud),
part_(part)
{}
// Member functions
inline label part() const
{
return part_;
}
};
private:
// Private data
//- Linear velocity of monoatomic
vector v_;
//- Total linear acceleration of monoatomic
vector a_;
//-
vector specialPosition_;
//-
scalar potentialEnergy_;
// - r_ij f_ij, stress dyad
tensor rf_;
// // - r_ij outer product f_ij: virial contribution
// tensor rDotf_;
//-
label special_;
//-
label id_;
//-
FixedList<vector, 1> siteForces_;
//-
FixedList<vector, 1> sitePositions_;
public:
//- Runtime type information
TypeName("monoatomic");
friend class Cloud<monoatomic>;
// Constructors
//- Construct with macroscopic description
inline monoatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const scalar temperature,
const vector& bulkVelocity,
const vector& specialPosition,
const constantProperties& constProps,
trackingData& td,
const label special,
const label id
);
//- Construct from all components
inline monoatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const vector& v,
const vector& a,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
const label id
);
//- Construct from Istream
monoatomic
(
const polyMesh& mesh,
Istream& is,
bool readFields = true
);
//- Construct and return a clone
autoPtr<particle> clone() const
{
return autoPtr<particle>(new monoatomic(*this));
}
//- Factory class to read-construct particles used for
// parallel transfer
class iNew
{
const polyMesh& mesh_;
public:
iNew(const polyMesh& mesh)
:
mesh_(mesh)
{}
autoPtr<monoatomic> operator()(Istream& is) const
{
return autoPtr<monoatomic>(new monoatomic(mesh_, is, true));
}
};
// Member Functions
// Tracking
//-
bool move(trackingData&, const scalar trackTime);
//-
virtual void transformProperties(const tensor& T);
//-
virtual void transformProperties(const vector& separation);
//-
void setSitePositions(const constantProperties& constProps);
//-
void setSiteSizes(label size);
// Access
//-
inline const vector& v() const;
//-
inline vector& v();
//-
inline const vector& a() const;
//-
inline vector& a();
//-
inline const FixedList<vector, 1>& siteForces() const;
//-
inline FixedList<vector, 1>& siteForces();
//-
inline const FixedList<vector, 1>& sitePositions() const;
//-
inline FixedList<vector, 1>& sitePositions();
//-
inline const vector& specialPosition() const;
//-
inline vector& specialPosition();
//-
inline scalar potentialEnergy() const;
//-
inline scalar& potentialEnergy();
//-
inline const tensor& rf() const;
//-
inline tensor& rf();
//-
inline label special() const;
//-
inline bool tethered() const;
//-
inline label id() const;
// Member Operators
//- Overridable function to handle the particle hitting a patch
// Executed before other patch-hitting functions
bool hitPatch
(
const polyPatch&,
trackingData& td,
const label patchI,
const scalar trackFraction,
const tetIndices& tetIs
);
//- Overridable function to handle the particle hitting a processorPatch
void hitProcessorPatch
(
const processorPolyPatch&,
trackingData& td
);
//- Overridable function to handle the particle hitting a wallPatch
void hitWallPatch
(
const wallPolyPatch&,
trackingData& td,
const tetIndices&
);
//- Overridable function to handle the particle hitting a polyPatch
void hitPatch
(
const polyPatch&,
trackingData& td
);
// I-O
//- Read
static void readFields(Cloud<monoatomic>& mC);
//- Write
static void writeFields(const Cloud<monoatomic>& mC);
//- Show info
static void info(trackingData& td);
// IOstream Operators
friend Ostream& operator<<(Ostream&, const monoatomic&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "monoatomicI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

328
src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicI.H
View File

@ -1,328 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::monoatomic::constantProperties::constantProperties()
:
sites_(),
electrostaticSites_(-1),
pairPotSites_(-1),
mass_(0)
{}
inline Foam::monoatomic::constantProperties::constantProperties
(
const dictionary& dict,
const List<label>& siteIds
)
:
sites_(),
electrostaticSites_(-1),
pairPotSites_(-1),
mass_(0)
{
if (siteIds.size() != 1)
{
FatalErrorIn
(
"Foam::monoatomic::constantProperties::constantProperties"
"("
"const dictionary& dict, "
"const List<label>& siteIds"
")"
)
<< "monoatomic, single site only, given: " << dict
<< nl << abort(FatalError);
}
FixedList<word, 1> siteIdNames = dict.lookup("siteIds");
FixedList<vector, 1> siteReferencePositions
(
dict.lookup("siteReferencePositions")
);
FixedList<scalar, 1> siteMasses(dict.lookup("siteMasses"));
FixedList<scalar, 1> siteCharges(dict.lookup("siteCharges"));
FixedList<word, 1> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
constPropSite site = sites_[0];
site = constPropSite
(
siteReferencePositions[0],
siteMasses[0],
siteCharges[0],
siteIds[0],
siteIdNames[0],
(findIndex(pairPotentialIds, siteIdNames[0]) != -1), // pair
(mag(siteCharges[0]) > VSMALL) // charge
);
mass_ = site.siteMass();
if (site.pairPotentialSite())
{
pairPotSites_[0] = 0;
}
if (site.electrostaticSite())
{
electrostaticSites_[0] = 0;
}
if (!site.pairPotentialSite() && !site.electrostaticSite())
{
WarningIn
(
"Foam::monoatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
)
<< siteIdNames[0] << " is a non-interacting site." << endl;
}
// Single site monoatomic - no rotational motion.
site.siteReferencePosition() = vector::zero;
}
inline Foam::monoatomic::monoatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const scalar temperature,
const vector& bulkVelocity,
const vector& specialPosition,
const constantProperties& cP,
monoatomic::trackingData& td,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
v_(bulkVelocity),
a_(vector::zero),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(vector::zero),
sitePositions_()
{
setSitePositions(cP);
v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
}
inline Foam::monoatomic::monoatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const vector& v,
const vector& a,
const vector& specialPosition,
const constantProperties& cP,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
v_(v),
a_(a),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(vector::zero),
sitePositions_()
{
setSitePositions(cP);
}
// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
inline const Foam::FixedList<Foam::constPropSite, 1>&
Foam::monoatomic::constantProperties::sites() const
{
return sites_;
}
inline const Foam::FixedList<Foam::label, 1>&
Foam::monoatomic::constantProperties::pairPotSites() const
{
return pairPotSites_;
}
inline const Foam::FixedList<Foam::label, 1>&
Foam::monoatomic::constantProperties::electrostaticSites() const
{
return electrostaticSites_;
}
inline Foam::label
Foam::monoatomic::constantProperties::degreesOfFreedom() const
{
return 3;
}
inline Foam::scalar Foam::monoatomic::constantProperties::mass() const
{
return mass_;
}
inline Foam::label Foam::monoatomic::constantProperties::nSites() const
{
return 1;
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline const Foam::vector& Foam::monoatomic::v() const
{
return v_;
}
inline Foam::vector& Foam::monoatomic::v()
{
return v_;
}
inline const Foam::vector& Foam::monoatomic::a() const
{
return a_;
}
inline Foam::vector& Foam::monoatomic::a()
{
return a_;
}
inline const Foam::FixedList<Foam::vector, 1>&
Foam::monoatomic::siteForces() const
{
return siteForces_;
}
inline Foam::FixedList<Foam::vector, 1>& Foam::monoatomic::siteForces()
{
return siteForces_;
}
inline const Foam::FixedList<Foam::vector, 1>&
Foam::monoatomic::sitePositions() const
{
return sitePositions_;
}
inline Foam::FixedList<Foam::vector, 1>& Foam::monoatomic::sitePositions()
{
return sitePositions_;
}
inline const Foam::vector& Foam::monoatomic::specialPosition() const
{
return specialPosition_;
}
inline Foam::vector& Foam::monoatomic::specialPosition()
{
return specialPosition_;
}
inline Foam::scalar Foam::monoatomic::potentialEnergy() const
{
return potentialEnergy_;
}
inline Foam::scalar& Foam::monoatomic::potentialEnergy()
{
return potentialEnergy_;
}
inline const Foam::tensor& Foam::monoatomic::rf() const
{
return rf_;
}
inline Foam::tensor& Foam::monoatomic::rf()
{
return rf_;
}
inline Foam::label Foam::monoatomic::special() const
{
return special_;
}
inline bool Foam::monoatomic::tethered() const
{
return special_ == SPECIAL_TETHERED;
}
inline Foam::label Foam::monoatomic::id() const
{
return id_;
}
// ************************************************************************* //

305
src/lagrangian/molecularDynamics/molecule/molecules/monoatomic/monoatomicIO.C
View File

@ -1,305 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "monoatomic.H"
#include "IOstreams.H"
#include "Cloud.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::monoatomic::monoatomic
(
const polyMesh& mesh,
Istream& is,
bool readFields
)
:
particle(mesh, is, readFields),
v_(vector::zero),
a_(vector::zero),
specialPosition_(vector::zero),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(0),
id_(0),
siteForces_(),
sitePositions_()
{
if (readFields)
{
if (is.format() == IOstream::ASCII)
{
is >> v_;
is >> a_;
is >> siteForces_;
potentialEnergy_ = readScalar(is);
is >> rf_;
special_ = readLabel(is);
id_ = readLabel(is);
is >> sitePositions_;
is >> specialPosition_;
}
else
{
is.read
(
reinterpret_cast<char*>(&v_),
sizeof(v_)
+ sizeof(a_)
+ sizeof(specialPosition_)
+ sizeof(potentialEnergy_)
+ sizeof(rf_)
+ sizeof(special_)
+ sizeof(id_)
);
is >> siteForces_ >> sitePositions_;
}
}
// Check state of Istream
is.check
(
"Foam::monoatomic::monoatomic"
"("
"const polyMesh& mesh,"
"Istream& is,"
"bool readFields"
")"
);
}
void Foam::monoatomic::readFields(Cloud<monoatomic>& mC)
{
if (!mC.size())
{
return;
}
particle::readFields(mC);
IOField<vector> v(mC.fieldIOobject("v", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, v);
IOField<vector> a(mC.fieldIOobject("a", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, a);
IOField<vector> specialPosition
(
mC.fieldIOobject("specialPosition", IOobject::MUST_READ)
);
mC.checkFieldIOobject(mC, specialPosition);
IOField<label> special(mC.fieldIOobject("special", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, special);
IOField<label> id(mC.fieldIOobject("id", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, id);
label i = 0;
forAllIter(typename Cloud<monoatomic>, mC, iter)
{
monoatomic& mol = iter();
mol.v_ = v[i];
mol.a_ = a[i];
mol.specialPosition_ = specialPosition[i];
mol.special_ = special[i];
mol.id_ = id[i];
i++;
}
}
void Foam::monoatomic::writeFields(const Cloud<monoatomic>& mC)
{
particle::writeFields(mC);
label np = mC.size();
IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
IOField<vector> a(mC.fieldIOobject("a", IOobject::NO_READ), np);
IOField<vector> specialPosition
(
mC.fieldIOobject("specialPosition", IOobject::NO_READ),
np
);
IOField<label> special(mC.fieldIOobject("special", IOobject::NO_READ), np);
IOField<label> id(mC.fieldIOobject("id", IOobject::NO_READ), np);
label i = 0;
forAllConstIter(typename Cloud<monoatomic>, mC, iter)
{
const monoatomic& mol = iter();
v[i] = mol.v_;
a[i] = mol.a_;
specialPosition[i] = mol.specialPosition_;
special[i] = mol.special_;
id[i] = mol.id_;
i++;
}
v.write();
a.write();
specialPosition.write();
special.write();
id.write();
}
void Foam::monoatomic::info(monoatomic::trackingData& td)
{
vector totalLinearMomentum(vector::zero);
scalar maxVelocityMag = 0.0;
scalar totalMass = 0.0;
scalar totalLinearKE = 0.0;
scalar totalPE = 0.0;
scalar totalrDotf = 0.0;
label nMols = td.cloud().size();
forAllConstIter(typename Cloud<monoatomic>, td.cloud(), mol)
{
const label molId = mol().id();
scalar molMass(td.cloud().constProps(molId).mass());
totalMass += molMass;
}
forAllConstIter(typename Cloud<monoatomic>, td.cloud(), mol)
{
const label molId = mol().id();
const monoatomic::constantProperties cP
(
td.cloud().constProps(molId)
);
scalar molMass(cP.mass());
const vector& molV(mol().v());
totalLinearMomentum += molV * molMass;
if (mag(molV) > maxVelocityMag)
{
maxVelocityMag = mag(molV);
}
totalLinearKE += 0.5*molMass*magSqr(molV);
totalPE += mol().potentialEnergy();
totalrDotf += tr(mol().rf());
}
scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
if (Pstream::parRun())
{
reduce(totalLinearMomentum, sumOp<vector>());
reduce(maxVelocityMag, maxOp<scalar>());
reduce(totalMass, sumOp<scalar>());
reduce(totalLinearKE, sumOp<scalar>());
reduce(totalPE, sumOp<scalar>());
reduce(totalrDotf, sumOp<scalar>());
reduce(nMols, sumOp<label>());
reduce(meshVolume, sumOp<scalar>());
}
if (nMols)
{
Info<< nl << "Number of molecules in " << td.cloud().name() << " = "
<< nMols << nl
<< " Overall number density = "
<< nMols/meshVolume << nl
<< " Overall mass density = "
<< totalMass/meshVolume << nl
<< " Average linear momentum per molecule = "
<< totalLinearMomentum/nMols << ' '
<< mag(totalLinearMomentum)/nMols << nl
<< " maximum |velocity| = "
<< maxVelocityMag << nl
<< " Average linear KE per molecule = "
<< totalLinearKE/nMols << nl
<< " Average angular KE per molecule = "
<< totalPE/nMols << nl
<< " Average TE per molecule = "
<<
(
totalLinearKE
+ totalPE
)
/nMols
<< nl << endl;
}
else
{
Info<< nl << "No molecules in " << td.cloud().name() << endl;
}
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const monoatomic& mol)
{
if (os.format() == IOstream::ASCII)
{
os << token::SPACE << static_cast<const particle&>(mol)
<< token::SPACE << mol.face()
<< token::SPACE << mol.stepFraction()
<< token::SPACE << mol.v_
<< token::SPACE << mol.a_
<< token::SPACE << mol.specialPosition_
<< token::SPACE << mol.potentialEnergy_
<< token::SPACE << mol.rf_
<< token::SPACE << mol.special_
<< token::SPACE << mol.id_
<< token::SPACE << mol.siteForces_
<< token::SPACE << mol.sitePositions_;
}
else
{
os << static_cast<const particle&>(mol);
os.write
(
reinterpret_cast<const char*>(&mol.v_),
sizeof(mol.v_)
+ sizeof(mol.a_)
+ sizeof(mol.specialPosition_)
+ sizeof(mol.potentialEnergy_)
+ sizeof(mol.rf_)
+ sizeof(mol.special_)
+ sizeof(mol.id_)
);
os << mol.siteForces_ << mol.sitePositions_;
}
// Check state of Ostream
os.check
(
"Foam::Ostream& Foam::operator<<"
"(Foam::Ostream&, const Foam::monoatomic&)"
);
return os;
}
// ************************************************************************* //

320
src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
View File

@ -1,320 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*----------------------------------------------------------------------------*/
#include "polyatomic.H"
#include "Random.H"
#include "Time.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(polyatomic, 0);
defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tensor Foam::polyatomic::rotationTensorX(scalar phi) const
{
return tensor
(
1, 0, 0,
0, Foam::cos(phi), -Foam::sin(phi),
0, Foam::sin(phi), Foam::cos(phi)
);
}
Foam::tensor Foam::polyatomic::rotationTensorY(scalar phi) const
{
return tensor
(
Foam::cos(phi), 0, Foam::sin(phi),
0, 1, 0,
-Foam::sin(phi), 0, Foam::cos(phi)
);
}
Foam::tensor Foam::polyatomic::rotationTensorZ(scalar phi) const
{
return tensor
(
Foam::cos(phi), -Foam::sin(phi), 0,
Foam::sin(phi), Foam::cos(phi), 0,
0, 0, 1
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::polyatomic::move
(
polyatomic::trackingData& td,
const scalar trackTime
)
{
td.switchProcessor = false;
td.keepParticle = true;
if (special_ != SPECIAL_FROZEN)
{
return td.keepParticle;
}
const constantProperties& constProps(td.cloud().constProps(id_));
if (td.part() == trackingData::tpFirstVelocityHalfStep)
{
// First leapfrog velocity adjust part, required before tracking+force
// part
v_ += 0.5*trackTime*a_;
pi_ += 0.5*trackTime*tau_;
}
else if (td.part() == trackingData::tpLinearTrack)
{
// Leapfrog tracking part
scalar tEnd = (1.0 - stepFraction())*trackTime;
scalar dtMax = tEnd;
while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
{
// set the lagrangian time-step
scalar dt = min(dtMax, tEnd);
dt *= trackToFace(position() + dt*v_, td);
tEnd -= dt;
stepFraction() = 1.0 - tEnd/trackTime;
}
}
else if (td.part() == trackingData::tpRotationalTrack)
{
// Leapfrog orientation adjustment, carried out before force calculation
// but after tracking stage, i.e. rotation carried once linear motion
// complete.
if (!constProps.pointMolecule())
{
const diagTensor& momentOfInertia(constProps.momentOfInertia());
tensor R;
if (!constProps.linearMolecule())
{
R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
pi_ = pi_ & R;
Q_ = Q_ & R;
}
R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
pi_ = pi_ & R;
Q_ = Q_ & R;
R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
pi_ = pi_ & R;
Q_ = Q_ & R;
R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
pi_ = pi_ & R;
Q_ = Q_ & R;
if (!constProps.linearMolecule())
{
R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
pi_ = pi_ & R;
Q_ = Q_ & R;
}
}
setSitePositions(constProps);
}
else if (td.part() == trackingData::tpSecondVelocityHalfStep)
{
// Second leapfrog velocity adjust part, required after tracking+force
// part
scalar m = constProps.mass();
a_ = vector::zero;
tau_ = vector::zero;
forAll(siteForces_, sI)
{
const vector& f = siteForces_[sI];
a_ += f/m;
tau_ +=
constProps.sites()[sI].siteReferencePosition()
^ (Q_.T() & f);
}
v_ += 0.5*trackTime*a_;
pi_ += 0.5*trackTime*tau_;
if (constProps.pointMolecule())
{
tau_ = vector::zero;
pi_ = vector::zero;
}
if (constProps.linearMolecule())
{
tau_.x() = 0.0;
pi_.x() = 0.0;
}
}
else
{
FatalErrorIn("polyatomic::move(trackingData&, const scalar)") << nl
<< td.part() << " is an invalid part of the integration method."
<< abort(FatalError);
}
return td.keepParticle;
}
void Foam::polyatomic::transformProperties(const tensor& T)
{
particle::transformProperties(T);
Q_ = T & Q_;
v_ = transform(T, v_);
a_ = transform(T, a_);
pi_ = Q_.T() & transform(T, Q_ & pi_);
tau_ = Q_.T() & transform(T, Q_ & tau_);
rf_ = transform(T, rf_);
sitePositions_ = position_ + (T & (sitePositions_ - position_));
siteForces_ = T & siteForces_;
}
void Foam::polyatomic::transformProperties(const vector& separation)
{
particle::transformProperties(separation);
if (special_ == SPECIAL_TETHERED)
{
specialPosition_ += separation;
}
sitePositions_ = sitePositions_ + separation;
}
void Foam::polyatomic::setSitePositions(const constantProperties& constProps)
{
forAll(constProps.sites(), sI)
{
sitePositions_[sI] =
position_
+ (Q_ & constProps.sites()[sI].siteReferencePosition());
}
}
void Foam::polyatomic::setSiteSizes(label size)
{
sitePositions_.setSize(size);
siteForces_.setSize(size);
}
bool Foam::polyatomic::hitPatch
(
const polyPatch&,
trackingData&,
const label,
const scalar,
const tetIndices&
)
{
return false;
}
void Foam::polyatomic::hitProcessorPatch
(
const processorPolyPatch&,
trackingData& td
)
{
td.switchProcessor = true;
}
void Foam::polyatomic::hitWallPatch
(
const wallPolyPatch& wpp,
trackingData& td,
const tetIndices& tetIs
)
{
// Use of the normal from tetIs is not required as
// hasWallImpactDistance is false.
vector nw = normal();
nw /= mag(nw);
scalar vn = v_ & nw;
// Specular reflection
if (vn > 0)
{
v_ -= 2*vn*nw;
}
}
void Foam::polyatomic::hitPatch
(
const polyPatch&,
trackingData& td
)
{
td.keepParticle = false;
}
// ************************************************************************* //

483
src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
View File

@ -1,483 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::polyatomic
Description
Foam::polyatomic
SourceFiles
polyatomicI.H
polyatomic.C
polyatomicIO.C
\*---------------------------------------------------------------------------*/
#ifndef polyatomic_H
#define polyatomic_H
#include "particle.H"
#include "IOstream.H"
#include "autoPtr.H"
#include "diagTensor.H"
#include "constPropSite.H"
#include "MoleculeCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyatomic Declaration
\*---------------------------------------------------------------------------*/
class polyatomic
:
public particle
{
public:
// Values of special that are less than zero are for built-in functionality.
// Values greater than zero are user specifiable/expandable (i.e. test
// special_ >= SPECIAL_USER)
enum specialTypes
{
SPECIAL_TETHERED = -1,
SPECIAL_FROZEN = -2,
NOT_SPECIAL = 0,
SPECIAL_USER = 1
};
//- Class to hold polyatomic constant properties
class constantProperties
{
// Private data
//- Sites of mass, charge or interaction
List<constPropSite> sites_;
//- Which sites require electrostatic interactions
List<label> electrostaticSites_;
//- Which sites require pair interactions
List<label> pairPotSites_;
//- Moment of intertia (in principal axis configiration)
diagTensor momentOfInertia_;
//-
scalar mass_;
// Private Member Functions
//-
bool linearMoleculeTest() const;
public:
//-
inline constantProperties();
//- Construct from dictionary
inline constantProperties
(
const dictionary& dict,
const labelList& siteIds
);
// Member functions
//-
inline const List<constPropSite>& sites() const;
//-
inline const List<label>& pairPotSites() const;
//-
inline const List<label>& electrostaticSites() const;
//-
inline const diagTensor& momentOfInertia() const;
//-
inline bool linearMolecule() const;
//-
inline bool pointMolecule() const;
//-
inline label degreesOfFreedom() const;
//-
inline scalar mass() const;
//-
inline label nSites() const;
};
//- Class used to pass tracking data to the trackToFace function
class trackingData
:
public particle::TrackingData<MoleculeCloud<polyatomic> >
{
public:
enum trackPart
{
tpFirstVelocityHalfStep,
tpLinearTrack,
tpRotationalTrack,
tpSecondVelocityHalfStep,
tpAccess
};
private:
// label specifying which part of the integration algorithm is taking
label part_;
public:
// Constructors
trackingData(MoleculeCloud<polyatomic>& cloud, trackPart part)
:
particle::TrackingData<MoleculeCloud<polyatomic> >(cloud),
part_(part)
{}
// Member functions
inline label part() const
{
return part_;
}
};
private:
// Private data
//- Orientation, stored as the rotation tensor to transform
// from the polyatomic to the global reference frame, i.e.:
// globalVector = Q_ & bodyLocalVector
// bodyLocalVector = Q_.T() & globalVector
tensor Q_;
//- Linear velocity of polyatomic
vector v_;
//- Total linear acceleration of polyatomic
vector a_;
//- Angular momentum of polyatomic, in body local reference frame
vector pi_;
//- Total torque on polyatomic, in body local reference frame
vector tau_;
//-
vector specialPosition_;
//-
scalar potentialEnergy_;
// - r_ij f_ij, stress dyad
tensor rf_;
// // - r_ij outer product f_ij: virial contribution
// tensor rDotf_;
//-
label special_;
//-
label id_;
//-
List<vector> siteForces_;
//-
List<vector> sitePositions_;
// Private Member Functions
//-
tensor rotationTensorX(scalar deltaT) const;
//-
tensor rotationTensorY(scalar deltaT) const;
//-
tensor rotationTensorZ(scalar deltaT) const;
public:
//- Runtime type information
TypeName("polyatomic");
friend class Cloud<polyatomic>;
// Constructors
//- Construct with macroscopic description
inline polyatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const scalar temperature,
const vector& bulkVelocity,
const vector& specialPosition,
const constantProperties& constProps,
polyatomic::trackingData& td,
const label special,
const label id
);
//- Construct from all components
inline polyatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const label special,
const label id
);
//- Construct from Istream
polyatomic
(
const polyMesh& mesh,
Istream& is,
bool readFields = true
);
//- Construct and return a clone
autoPtr<particle> clone() const
{
return autoPtr<particle>(new polyatomic(*this));
}
//- Factory class to read-construct particles used for
// parallel transfer
class iNew
{
const polyMesh& mesh_;
public:
iNew(const polyMesh& mesh)
:
mesh_(mesh)
{}
autoPtr<polyatomic> operator()(Istream& is) const
{
return autoPtr<polyatomic>(new polyatomic(mesh_, is, true));
}
};
// Member Functions
// Tracking
//-
bool move(trackingData&, const scalar trackTime);
//-
virtual void transformProperties(const tensor& T);
//-
virtual void transformProperties(const vector& separation);
//-
void setSitePositions(const constantProperties& constProps);
//-
void setSiteSizes(label size);
// Access
//-
inline const tensor& Q() const;
//-
inline tensor& Q();
//-
inline const vector& v() const;
//-
inline vector& v();
//-
inline const vector& a() const;
//-
inline vector& a();
//-
inline const vector& pi() const;
//-
inline vector& pi();
//-
inline const vector& tau() const;
//-
inline vector& tau();
//-
inline const List<vector>& siteForces() const;
//-
inline List<vector>& siteForces();
//-
inline const List<vector>& sitePositions() const;
//-
inline List<vector>& sitePositions();
//-
inline const vector& specialPosition() const;
//-
inline vector& specialPosition();
//-
inline scalar potentialEnergy() const;
//-
inline scalar& potentialEnergy();
//-
inline const tensor& rf() const;
//-
inline tensor& rf();
//-
inline label special() const;
//-
inline bool tethered() const;
//-
inline label id() const;
// Member Operators
//- Overridable function to handle the particle hitting a patch
// Executed before other patch-hitting functions
bool hitPatch
(
const polyPatch&,
trackingData& td,
const label patchI,
const scalar trackFraction,
const tetIndices& tetIs
);
//- Overridable function to handle the particle hitting a processorPatch
void hitProcessorPatch
(
const processorPolyPatch&,
trackingData& td
);
//- Overridable function to handle the particle hitting a wallPatch
void hitWallPatch
(
const wallPolyPatch&,
trackingData& td,
const tetIndices&
);
//- Overridable function to handle the particle hitting a polyPatch
void hitPatch
(
const polyPatch&,
trackingData& td
);
// I-O
//- Read
static void readFields(Cloud<polyatomic>& mC);
//- Write
static void writeFields(const Cloud<polyatomic>& mC);
//- Show info
static void info(trackingData& td);
// IOstream Operators
friend Ostream& operator<<(Ostream&, const polyatomic&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "polyatomicI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

653
src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicI.H
View File

@ -1,653 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::polyatomic::constantProperties::constantProperties()
:
sites_(),
electrostaticSites_(),
pairPotSites_(),
momentOfInertia_(diagTensor(0, 0, 0)),
mass_(0)
{}
inline Foam::polyatomic::constantProperties::constantProperties
(
const dictionary& dict,
const labelList& siteIds
)
:
sites_(),
electrostaticSites_(),
pairPotSites_(),
momentOfInertia_(),
mass_(0)
{
List<word> siteIdNames = dict.lookup("siteIds");
List<vector> siteReferencePositions(dict.lookup("siteReferencePositions"));
List<scalar> siteMasses(dict.lookup("siteMasses"));
List<scalar> siteCharges(dict.lookup("siteCharges"));
List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
sites_.setSize(siteReferencePositions.size());
if
(
(siteIdNames.size() != sites_.size())
|| (siteCharges.size() != sites_.size())
)
{
FatalErrorIn
(
"Foam::polyatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
)
<< "Sizes of site id, charge and "
<< "referencePositions are not the same: " << nl
<< siteIdNames << nl
<< siteReferencePositions << nl
<< siteCharges << nl
<< abort(FatalError);
}
electrostaticSites_.setSize(sites_.size(), -1);
pairPotSites_.setSize(sites_.size(), -1);
label pairI = 0;
label electrostaticI = 0;
forAll(sites_, sI)
{
sites_[sI] = constPropSite
(
siteReferencePositions[sI],
siteMasses[sI],
siteCharges[sI],
siteIds[sI],
siteIdNames[sI],
(findIndex(pairPotentialIds, siteIdNames[sI]) != -1), // pair
(mag(siteCharges[sI]) > VSMALL) // charge
);
if (sites_[sI].pairPotentialSite())
{
pairPotSites_[pairI++] = sI;
}
if (sites_[sI].electrostaticSite())
{
electrostaticSites_[electrostaticI++] = sI;
}
if (!sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
{
WarningIn
(
"Foam::polyatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
)
<< siteIdNames[sI] << " is a non-interacting site." << endl;
}
}
pairPotSites_.setSize(pairI);
electrostaticSites_.setSize(electrostaticI);
// Calculate the centre of mass of the body and subtract it from each
// position
vector centreOfMass(vector::zero);
forAll(sites_, sI)
{
mass_ += sites_[sI].siteMass();
centreOfMass +=
sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
}
centreOfMass /= mass_;
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() -= centreOfMass;
}
if (sites_.size() == 1)
{
// Single site polyatomic - no rotational motion.
sites_[0].siteReferencePosition() = vector::zero;
momentOfInertia_ = diagTensor(-1, -1, -1);
}
else if (linearMoleculeTest())
{
// Linear polyatomic.
Info<< nl << "Linear polyatomic." << endl;
vector dir =
sites_[1].siteReferencePosition()
- sites_[0].siteReferencePosition();
dir /= mag(dir);
tensor Q = rotationTensor(dir, vector(1,0,0));
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() =
(Q & sites_[sI].siteReferencePosition());
}
// The rotation was around the centre of mass but remove any
// components that have crept in due to floating point errors
centreOfMass = vector::zero;
forAll(sites_, sI)
{
centreOfMass +=
sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
}
centreOfMass /= mass_;
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() -= centreOfMass;
}
diagTensor momOfInertia = diagTensor::zero;
forAll(sites_, sI)
{
const vector& p(sites_[sI].siteReferencePosition());
momOfInertia +=
sites_[sI].siteMass()*diagTensor(0, p.x()*p.x(), p.x()*p.x());
}
momentOfInertia_ = diagTensor
(
-1,
momOfInertia.yy(),
momOfInertia.zz()
);
}
else
{
// Fully 6DOF polyatomic
// Calculate the inertia tensor in the current orientation
tensor momOfInertia(tensor::zero);
forAll(sites_, sI)
{
const vector& p(sites_[sI].siteReferencePosition());
momOfInertia += sites_[sI].siteMass()*tensor
(
p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
);
}
if (eigenValues(momOfInertia).x() < VSMALL)
{
FatalErrorIn("polyatomic::constantProperties::constantProperties")
<< "An eigenvalue of the inertia tensor is zero, the "
<< "polyatomic " << dict.name() << " is not a valid 6DOF shape."
<< nl << abort(FatalError);
}
// Normalise the inertia tensor magnitude to avoid SMALL numbers in the
// components causing problems
momOfInertia /= eigenValues(momOfInertia).x();
tensor e = eigenVectors(momOfInertia);
// Calculate the transformation between the principle axes to the
// global axes
tensor Q =
vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
// Transform the site positions
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() =
(Q & sites_[sI].siteReferencePosition());
}
// Recalculating the moment of inertia with the new site positions
// The rotation was around the centre of mass but remove any
// components that have crept in due to floating point errors
centreOfMass = vector::zero;
forAll(sites_, sI)
{
centreOfMass +=
sites_[sI].siteReferencePosition()*sites_[sI].siteMass();
}
centreOfMass /= mass_;
forAll(sites_, sI)
{
sites_[sI].siteReferencePosition() -= centreOfMass;
}
// Calculate the moment of inertia in the principle component
// reference frame
momOfInertia = tensor::zero;
forAll(sites_, sI)
{
const vector& p(sites_[sI].siteReferencePosition());
momOfInertia += sites_[sI].siteMass()*tensor
(
p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
-p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
-p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y()
);
}
momentOfInertia_ = diagTensor
(
momOfInertia.xx(),
momOfInertia.yy(),
momOfInertia.zz()
);
}
}
inline Foam::polyatomic::polyatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const scalar temperature,
const vector& bulkVelocity,
const vector& specialPosition,
const constantProperties& cP,
polyatomic::trackingData& td,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
Q_(I),
v_(bulkVelocity),
a_(vector::zero),
pi_(vector::zero),
tau_(vector::zero),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(cP.nSites(), vector::zero),
sitePositions_(cP.nSites())
{
setSitePositions(cP);
v_ += td.cloud().equipartitionLinearVelocity(temperature, cP.mass());
if (!cP.pointMolecule())
{
pi_ = td.cloud().equipartitionAngularMomentum(temperature, cP);
scalar phi(td.cloud().rndGen().scalar01()*twoPi);
scalar theta(td.cloud().rndGen().scalar01()*twoPi);
scalar psi(td.cloud().rndGen().scalar01()*twoPi);
Q_ = tensor
(
cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
sin(psi)*sin(theta),
- sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
- sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
cos(psi)*sin(theta),
sin(theta)*sin(phi),
- sin(theta)*cos(phi),
cos(theta)
);
}
}
Foam::polyatomic::polyatomic
(
const polyMesh& mesh,
const vector& position,
const label cellI,
const label tetFaceI,
const label tetPtI,
const tensor& Q,
const vector& v,
const vector& a,
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& cP,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
Q_(Q),
v_(v),
a_(a),
pi_(pi),
tau_(tau),
specialPosition_(specialPosition),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(cP.nSites(), vector::zero),
sitePositions_(cP.nSites())
{
setSitePositions(cP);
}
// * * * * * * * constantProperties Private Member Functions * * * * * * * * //
inline bool Foam::polyatomic::constantProperties::linearMoleculeTest() const
{
if (sites_.size() == 2)
{
return true;
}
vector refDir =
sites_[1].siteReferencePosition()
- sites_[0].siteReferencePosition();
refDir /= mag(refDir);
for
(
label i = 2;
i < sites_.size();
i++
)
{
vector dir =
sites_[i].siteReferencePosition()
- sites_[i - 1].siteReferencePosition();
dir /= mag(dir);
if (mag(refDir & dir) < 1 - SMALL)
{
return false;
}
}
return true;
}
// * * * * * * * constantProperties Member Functions * * * * * * * * * * * * //
inline const Foam::List<Foam::constPropSite>&
Foam::polyatomic::constantProperties::sites() const
{
return sites_;
}
inline const Foam::List<Foam::label>&
Foam::polyatomic::constantProperties::pairPotSites() const
{
return pairPotSites_;
}
inline const Foam::List<Foam::label>&
Foam::polyatomic::constantProperties::electrostaticSites() const
{
return electrostaticSites_;
}
inline const Foam::diagTensor&
Foam::polyatomic::constantProperties::momentOfInertia() const
{
return momentOfInertia_;
}
inline bool Foam::polyatomic::constantProperties::linearMolecule() const
{
return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
}
inline bool Foam::polyatomic::constantProperties::pointMolecule() const
{
return (momentOfInertia_.zz() < 0);
}
inline Foam::label
Foam::polyatomic::constantProperties::degreesOfFreedom() const
{
if (linearMolecule())
{
return 5;
}
else if (pointMolecule())
{
return 3;
}
else
{
return 6;
}
}
inline Foam::scalar Foam::polyatomic::constantProperties::mass() const
{
return mass_;
}
inline Foam::label Foam::polyatomic::constantProperties::nSites() const
{
return sites_.size();
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline const Foam::tensor& Foam::polyatomic::Q() const
{
return Q_;
}
inline Foam::tensor& Foam::polyatomic::Q()
{
return Q_;
}
inline const Foam::vector& Foam::polyatomic::v() const
{
return v_;
}
inline Foam::vector& Foam::polyatomic::v()
{
return v_;
}
inline const Foam::vector& Foam::polyatomic::a() const
{
return a_;
}
inline Foam::vector& Foam::polyatomic::a()
{
return a_;
}
inline const Foam::vector& Foam::polyatomic::pi() const
{
return pi_;
}
inline Foam::vector& Foam::polyatomic::pi()
{
return pi_;
}
inline const Foam::vector& Foam::polyatomic::tau() const
{
return tau_;
}
inline Foam::vector& Foam::polyatomic::tau()
{
return tau_;
}
inline const Foam::List<Foam::vector>& Foam::polyatomic::siteForces() const
{
return siteForces_;
}
inline Foam::List<Foam::vector>& Foam::polyatomic::siteForces()
{
return siteForces_;
}
inline const Foam::List<Foam::vector>& Foam::polyatomic::sitePositions() const
{
return sitePositions_;
}
inline Foam::List<Foam::vector>& Foam::polyatomic::sitePositions()
{
return sitePositions_;
}
inline const Foam::vector& Foam::polyatomic::specialPosition() const
{
return specialPosition_;
}
inline Foam::vector& Foam::polyatomic::specialPosition()
{
return specialPosition_;
}
inline Foam::scalar Foam::polyatomic::potentialEnergy() const
{
return potentialEnergy_;
}
inline Foam::scalar& Foam::polyatomic::potentialEnergy()
{
return potentialEnergy_;
}
inline const Foam::tensor& Foam::polyatomic::rf() const
{
return rf_;
}
inline Foam::tensor& Foam::polyatomic::rf()
{
return rf_;
}
inline Foam::label Foam::polyatomic::special() const
{
return special_;
}
inline bool Foam::polyatomic::tethered() const
{
return special_ == SPECIAL_TETHERED;
}
inline Foam::label Foam::polyatomic::id() const
{
return id_;
}
// ************************************************************************* //

414
src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomicIO.C
View File

@ -1,414 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "polyatomic.H"
#include "IOstreams.H"
#include "Cloud.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyatomic::polyatomic
(
const polyMesh& mesh,
Istream& is,
bool readFields
)
:
particle(mesh, is, readFields),
Q_(tensor::zero),
v_(vector::zero),
a_(vector::zero),
pi_(vector::zero),
tau_(vector::zero),
specialPosition_(vector::zero),
potentialEnergy_(0.0),
rf_(tensor::zero),
special_(0),
id_(0),
siteForces_(0),
sitePositions_(0)
{
if (readFields)
{
if (is.format() == IOstream::ASCII)
{
is >> Q_;
is >> v_;
is >> a_;
is >> pi_;
is >> tau_;
is >> siteForces_;
potentialEnergy_ = readScalar(is);
is >> rf_;
special_ = readLabel(is);
id_ = readLabel(is);
is >> sitePositions_;
is >> specialPosition_;
}
else
{
is.read
(
reinterpret_cast<char*>(&Q_),
sizeof(Q_)
+ sizeof(v_)
+ sizeof(a_)
+ sizeof(pi_)
+ sizeof(tau_)
+ sizeof(specialPosition_)
+ sizeof(potentialEnergy_)
+ sizeof(rf_)
+ sizeof(special_)
+ sizeof(id_)
);
is >> siteForces_ >> sitePositions_;
}
}
// Check state of Istream
is.check
(
"Foam::polyatomic::polyatomic"
"("
"const polyMesh& mesh,"
"Istream& is,"
"bool readFields"
")"
);
}
void Foam::polyatomic::readFields(Cloud<polyatomic>& mC)
{
if (!mC.size())
{
return;
}
particle::readFields(mC);
IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, Q);
IOField<vector> v(mC.fieldIOobject("v", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, v);
IOField<vector> a(mC.fieldIOobject("a", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, a);
IOField<vector> pi(mC.fieldIOobject("pi", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, pi);
IOField<vector> tau(mC.fieldIOobject("tau", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, tau);
IOField<vector> specialPosition
(
mC.fieldIOobject("specialPosition", IOobject::MUST_READ)
);
mC.checkFieldIOobject(mC, specialPosition);
IOField<label> special(mC.fieldIOobject("special", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, special);
IOField<label> id(mC.fieldIOobject("id", IOobject::MUST_READ));
mC.checkFieldIOobject(mC, id);
label i = 0;
forAllIter(typename Cloud<polyatomic>, mC, iter)
{
polyatomic& mol = iter();
mol.Q_ = Q[i];
mol.v_ = v[i];
mol.a_ = a[i];
mol.pi_ = pi[i];
mol.tau_ = tau[i];
mol.specialPosition_ = specialPosition[i];
mol.special_ = special[i];
mol.id_ = id[i];
i++;
}
}
void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
{
particle::writeFields(mC);
label np = mC.size();
IOField<tensor> Q(mC.fieldIOobject("Q", IOobject::NO_READ), np);
IOField<vector> v(mC.fieldIOobject("v", IOobject::NO_READ), np);
IOField<vector> a(mC.fieldIOobject("a", IOobject::NO_READ), np);
IOField<vector> pi(mC.fieldIOobject("pi", IOobject::NO_READ), np);
IOField<vector> tau(mC.fieldIOobject("tau", IOobject::NO_READ), np);
IOField<vector> specialPosition
(
mC.fieldIOobject("specialPosition", IOobject::NO_READ),
np
);
IOField<label> special(mC.fieldIOobject("special", IOobject::NO_READ), np);
IOField<label> id(mC.fieldIOobject("id", IOobject::NO_READ), np);
// Post processing fields
IOField<vector> piGlobal
(
mC.fieldIOobject("piGlobal", IOobject::NO_READ),
np
);
IOField<vector> tauGlobal
(
mC.fieldIOobject("tauGlobal", IOobject::NO_READ),
np
);
IOField<vector> orientation1
(
mC.fieldIOobject("orientation1", IOobject::NO_READ),
np
);
IOField<vector> orientation2
(
mC.fieldIOobject("orientation2", IOobject::NO_READ),
np
);
IOField<vector> orientation3
(
mC.fieldIOobject("orientation3", IOobject::NO_READ),
np
);
label i = 0;
forAllConstIter(typename Cloud<polyatomic>, mC, iter)
{
const polyatomic& mol = iter();
Q[i] = mol.Q_;
v[i] = mol.v_;
a[i] = mol.a_;
pi[i] = mol.pi_;
tau[i] = mol.tau_;
specialPosition[i] = mol.specialPosition_;
special[i] = mol.special_;
id[i] = mol.id_;
piGlobal[i] = mol.Q_ & mol.pi_;
tauGlobal[i] = mol.Q_ & mol.tau_;
orientation1[i] = mol.Q_ & vector(1,0,0);
orientation2[i] = mol.Q_ & vector(0,1,0);
orientation3[i] = mol.Q_ & vector(0,0,1);
i++;
}
Q.write();
v.write();
a.write();
pi.write();
tau.write();
specialPosition.write();
special.write();
id.write();
piGlobal.write();
tauGlobal.write();
orientation1.write();
orientation2.write();
orientation3.write();
}
void Foam::polyatomic::info(polyatomic::trackingData& td)
{
vector totalLinearMomentum(vector::zero);
vector totalAngularMomentum(vector::zero);
scalar maxVelocityMag = 0.0;
scalar totalMass = 0.0;
scalar totalLinearKE = 0.0;
scalar totalAngularKE = 0.0;
scalar totalPE = 0.0;
scalar totalrDotf = 0.0;
//vector CentreOfMass(vector::zero);
label nMols = td.cloud().size();
label dofs = 0;
forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
{
const label molId = mol().id();
scalar molMass(td.cloud().constProps(molId).mass());
totalMass += molMass;
//CentreOfMass += mol().position()*molMass;
}
// if (nMols)
// {
// CentreOfMass /= totalMass;
// }
forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
{
const label molId = mol().id();
const polyatomic::constantProperties cP
(
td.cloud().constProps(molId)
);
scalar molMass(cP.mass());
const diagTensor& molMoI(cP.momentOfInertia());
const vector& molV(mol().v());
const vector& molOmega(inv(molMoI) & mol().pi());
vector molPiGlobal = mol().Q() & mol().pi();
totalLinearMomentum += molV * molMass;
totalAngularMomentum += molPiGlobal;
//+((mol().position() - CentreOfMass) ^ (molV * molMass));
if (mag(molV) > maxVelocityMag)
{
maxVelocityMag = mag(molV);
}
totalLinearKE += 0.5*molMass*magSqr(molV);
totalAngularKE += 0.5*(molOmega & molMoI & molOmega);
totalPE += mol().potentialEnergy();
totalrDotf += tr(mol().rf());
dofs += cP.degreesOfFreedom();
}
scalar meshVolume = sum(td.cloud().mesh().cellVolumes());
if (Pstream::parRun())
{
reduce(totalLinearMomentum, sumOp<vector>());
reduce(totalAngularMomentum, sumOp<vector>());
reduce(maxVelocityMag, maxOp<scalar>());
reduce(totalMass, sumOp<scalar>());
reduce(totalLinearKE, sumOp<scalar>());
reduce(totalAngularKE, sumOp<scalar>());
reduce(totalPE, sumOp<scalar>());
reduce(totalrDotf, sumOp<scalar>());
reduce(nMols, sumOp<label>());
reduce(dofs, sumOp<label>());
reduce(meshVolume, sumOp<scalar>());
}
if (nMols)
{
Info<< nl << "Number of molecules in " << td.cloud().name() << " = "
<< nMols << nl
<< " Overall number density = "
<< nMols/meshVolume << nl
<< " Overall mass density = "
<< totalMass/meshVolume << nl
<< " Average linear momentum per molecule = "
<< totalLinearMomentum/nMols << ' '
<< mag(totalLinearMomentum)/nMols << nl
<< " Average angular momentum per molecule = "
<< totalAngularMomentum << ' '
<< mag(totalAngularMomentum)/nMols << nl
<< " maximum |velocity| = "
<< maxVelocityMag << nl
<< " Average linear KE per molecule = "
<< totalLinearKE/nMols << nl
<< " Average angular KE per molecule = "
<< totalAngularKE/nMols << nl
<< " Average PE per molecule = "
<< totalPE/nMols << nl
<< " Average TE per molecule = "
<<
(
totalLinearKE
+ totalAngularKE
+ totalPE
)
/nMols
<< nl << endl;
}
else
{
Info<< nl << "No molecules in " << td.cloud().name() << endl;
}
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const polyatomic& mol)
{
if (os.format() == IOstream::ASCII)
{
os << token::SPACE << static_cast<const particle&>(mol)
<< token::SPACE << mol.face()
<< token::SPACE << mol.stepFraction()
<< token::SPACE << mol.Q_
<< token::SPACE << mol.v_
<< token::SPACE << mol.a_
<< token::SPACE << mol.pi_
<< token::SPACE << mol.tau_
<< token::SPACE << mol.specialPosition_
<< token::SPACE << mol.potentialEnergy_
<< token::SPACE << mol.rf_
<< token::SPACE << mol.special_
<< token::SPACE << mol.id_
<< token::SPACE << mol.siteForces_
<< token::SPACE << mol.sitePositions_;
}
else
{
os << static_cast<const particle&>(mol);
os.write
(
reinterpret_cast<const char*>(&mol.Q_),
sizeof(mol.Q_)
+ sizeof(mol.v_)
+ sizeof(mol.a_)
+ sizeof(mol.pi_)
+ sizeof(mol.tau_)
+ sizeof(mol.specialPosition_)
+ sizeof(mol.potentialEnergy_)
+ sizeof(mol.rf_)
+ sizeof(mol.special_)
+ sizeof(mol.id_)
);
os << mol.siteForces_ << mol.sitePositions_;
}
// Check state of Ostream
os.check
(
"Foam::Ostream& Foam::operator<<"
"(Foam::Ostream&, const Foam::polyatomic&)"
);
return os;
}
// ************************************************************************* //

6
src/turbulenceModels/derivedFvPatchFields/porousBafflePressure/porousBafflePressureFvPatchFields.C
View File

@ -66,6 +66,8 @@ void Foam::porousBafflePressureFvPatchField<Foam::scalar>::updateCoeffs()
scalarField Un(phip/patch().magSf()); scalarField Un(phip/patch().magSf());
scalarField magUn(mag(Un));
if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0)) if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
{ {
const incompressible::turbulenceModel& model = const incompressible::turbulenceModel& model =
@ -76,7 +78,7 @@ void Foam::porousBafflePressureFvPatchField<Foam::scalar>::updateCoeffs()
const scalarField nuEffw = model.nuEff()().boundaryField()[patchI]; const scalarField nuEffw = model.nuEff()().boundaryField()[patchI];
jump_ = -(I_*nuEffw*mag(Un) + D_*0.5*magSqr(Un)*length_); jump_ = -sign(Un)*(I_*nuEffw + D_*0.5*magUn*length_)*magUn;
} }
else else
{ {
@ -93,7 +95,7 @@ void Foam::porousBafflePressureFvPatchField<Foam::scalar>::updateCoeffs()
Un /= rhow; Un /= rhow;
jump_ = -(I_*muEffw*mag(Un) + D_*0.5*rhow*magSqr(Un)*length_); jump_ = -sign(Un)*(I_*muEffw + D_*0.5*rhow*magUn*length_)*magUn;
} }
if (debug) if (debug)