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 'cvm' of /home/noisy3/OpenFOAM/OpenFOAM-dev into cvm

Browse Source
This commit is contained in:
mattijs
2011-07-06 11:21:26 +01:00
parent d6eeb3ab89 3598a450fb
commit bbe91bb651
65 changed files with 5148 additions and 615 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
View File

@ -1,4 +1,7 @@
EXE_DEBUG = -DFULLDEBUG -g -O0
EXE_INC = \
${EXE_DEBUG} \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \

42
applications/solvers/discreteMethods/molecularDynamics/mdFoam/createFields.H
View File

@ -12,6 +12,44 @@
mesh
);
potential pot(mesh);
word polyatomicCloudName("polyatomicCloud");
polyatomicCloud molecules(mesh, pot);
potential polyPot
(
mesh,
IOdictionary
(
IOobject
(
polyatomicCloudName + "Properties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
polyatomicCloud polyatomics(polyatomicCloudName, mesh, polyPot);
word monoatomicCloudName("monoatomicCloud");
potential monoPot
(
mesh,
IOdictionary
(
IOobject
(
monoatomicCloudName + "Properties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
monoatomicCloud monoatomics(monoatomicCloudName, mesh, monoPot);

9
applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C
View File

@ -30,6 +30,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "monoatomicCloud.H"
#include "polyatomicCloud.H"
int main(int argc, char *argv[])
@ -47,9 +48,13 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << endl;
molecules.evolve();
monoatomics.evolve();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
polyatomics.evolve();
runTime.write();
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

0
applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/files → applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/files
View File

0
applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options → applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/Make/options
View File

0
applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C → applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/mdEquilibrationFoam.C
View File

0
applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H → applications/solvers/discreteMethods/molecularDynamics/old/mdEquilibrationFoam/readmdEquilibrationDict.H
View File

6
applications/utilities/mesh/manipulation/topoSet/topoSet.C
View File

@ -208,7 +208,8 @@ int main(int argc, char *argv[])
IOobject::MUST_READ
);
Info<< "Read set " << setName << " with size "
<< currentSet().size() << endl;
<< returnReduce(currentSet().size(), sumOp<label>())
<< endl;
}
@ -292,7 +293,8 @@ int main(int argc, char *argv[])
if (currentSet.valid())
{
Info<< " Set " << currentSet().name()
<< " now size " << currentSet().size()
<< " now size "
<< returnReduce(currentSet().size(), sumOp<label>())
<< endl;
}
}

4
applications/utilities/preProcessing/mdInitialise/Make/options
View File

@ -3,7 +3,6 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
@ -13,6 +12,5 @@ EXE_LIBS = \
-lfiniteVolume \
-llagrangian \
-lmolecule \
-lpotential \
-lmolecularMeasurements
-lpotential

99
applications/utilities/preProcessing/mdInitialise/mdInitialise.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,8 +26,9 @@ Description
\*---------------------------------------------------------------------------*/
#include "md.H"
#include "fvCFD.H"
#include "polyatomicCloud.H"
#include "monoatomicCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -50,11 +51,15 @@ int main(int argc, char *argv[])
)
);
IOdictionary idListDict
word polyCloudName("polyatomicCloud");
const dictionary& polyDict(mdInitialiseDict.subDict(polyCloudName));
IOdictionary polyIdListDict
(
IOobject
(
"idList",
polyCloudName + "_idList",
mesh.time().constant(),
mesh,
IOobject::NO_READ,
@ -62,26 +67,88 @@ int main(int argc, char *argv[])
)
);
potential pot(mesh, mdInitialiseDict, idListDict);
potential polyPot
(
mesh,
polyDict,
IOdictionary
(
IOobject
(
polyCloudName + "Properties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
),
polyIdListDict
);
moleculeCloud molecules(mesh, pot, mdInitialiseDict);
polyatomicCloud poly
(
polyCloudName,
mesh,
polyPot,
polyDict
);
label totalMolecules = molecules.size();
if (Pstream::parRun())
{
reduce(totalMolecules, sumOp<label>());
}
Info<< nl << "Total number of molecules added: " << totalMolecules
Info<< nl << returnReduce(poly.size(), sumOp<label>()) << " added to "
<< poly.name()
<< nl << endl;
IOstream::defaultPrecision(15);
word monoCloudName("monoatomicCloud");
const dictionary& monoDict(mdInitialiseDict.subDict(monoCloudName));
IOdictionary monoIdListDict
(
IOobject
(
monoCloudName + "_idList",
mesh.time().constant(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
)
);
potential monoPot
(
mesh,
monoDict,
IOdictionary
(
IOobject
(
monoCloudName + "Properties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
),
monoIdListDict
);
monoatomicCloud mono
(
monoCloudName,
mesh,
monoPot,
monoDict
);
Info<< nl << returnReduce(mono.size(), sumOp<label>()) << " added to "
<< mono.name()
<< nl << endl;
if (!mesh.write())
{
FatalErrorIn(args.executable())
<< "Failed writing moleculeCloud."
<< "Failed writing."
<< nl << exit(FatalError);
}

38
src/OSspecific/POSIX/POSIX.C
View File

@ -57,6 +57,10 @@ Description
#include <netinet/in.h>
#ifdef USE_RANDOM
# include <climits>
#endif
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::POSIX, 0);
@ -68,16 +72,19 @@ pid_t Foam::pid()
return ::getpid();
}
pid_t Foam::ppid()
{
return ::getppid();
}
pid_t Foam::pgid()
{
return ::getpgrp();
}
bool Foam::env(const word& envName)
{
return ::getenv(envName.c_str()) != NULL;
@ -890,7 +897,6 @@ bool Foam::mvBak(const fileName& src, const std::string& ext)
}
// Remove a file, returning true if successful otherwise false
bool Foam::rm(const fileName& file)
{
@ -1221,4 +1227,34 @@ Foam::fileNameList Foam::dlLoaded()
}
void Foam::osRandomSeed(const label seed)
{
#ifdef USE_RANDOM
srandom((unsigned int)seed);
#else
srand48(seed);
#endif
}
Foam::label Foam::osRandomInteger()
{
#ifdef USE_RANDOM
return random();
#else
return lrand48();
#endif
}
Foam::scalar Foam::osRandomDouble()
{
#ifdef USE_RANDOM
return (scalar)random();
#else
return drand48();
#endif
}
// ************************************************************************* //

12
src/OpenFOAM/include/OSspecific.H
View File

@ -200,6 +200,18 @@ bool dlSymFound(void* handle, const std::string& symbol);
fileNameList dlLoaded();
// Low level random numbers. Use Random class instead.
//- Seed random number generator.
void osRandomSeed(const label seed);
//- Return random integer (uniform distribution between 0 and 2^31)
label osRandomInteger();
//- Return random double precision (uniform distribution between 0 and 1)
scalar osRandomDouble();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

33
src/OpenFOAM/primitives/random/Random/Random.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "Random.H"
#include "OSspecific.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -37,11 +38,6 @@ namespace Foam
# error "The random number generator may not work!"
#endif
#ifdef USE_RANDOM
# include <climits>
#else
# include <cstdlib>
#endif
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -57,22 +53,13 @@ Random::Random(const label seed)
Seed = 1;
}
# ifdef USE_RANDOM
srandom((unsigned int)Seed);
# else
srand48(Seed);
# endif
osRandomSeed(Seed);
}
int Random::bit()
{
# ifdef USE_RANDOM
if (random() > INT_MAX/2)
# else
if (lrand48() > INT_MAX/2)
# endif
if (osRandomInteger() > INT_MAX/2)
{
return 1;
}
@ -85,11 +72,7 @@ int Random::bit()
scalar Random::scalar01()
{
# ifdef USE_RANDOM
return (scalar)random()/INT_MAX;
# else
return drand48();
# endif
return osRandomDouble();
}
@ -140,11 +123,7 @@ tensor Random::tensor01()
label Random::integer(const label lower, const label upper)
{
# ifdef USE_RANDOM
return lower + (random() % (upper+1-lower));
# else
return lower + (lrand48() % (upper+1-lower));
# endif
return lower + (osRandomInteger() % (upper+1-lower));
}

10
src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,7 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "cachedRandom.H"
#include <cstdlib>
#include "OSspecific.H"
#if INT_MAX != 2147483647
# error "INT_MAX != 2147483647"
@ -37,7 +37,7 @@ Foam::scalar Foam::cachedRandom::scalar01()
{
if (sampleI_ < 0)
{
return drand48();
return osRandomDouble();
}
if (sampleI_ == samples_.size() - 1)
@ -76,7 +76,7 @@ Foam::cachedRandom::cachedRandom(const label seed, const label count)
}
// Initialise samples
srand48(seed_);
osRandomSeed(seed_);
forAll(samples_, i)
{
samples_[i] = drand48();
@ -98,7 +98,7 @@ Foam::cachedRandom::cachedRandom(const cachedRandom& cr, const bool reset)
) << "Copy constructor called, but samples not being cached. "
<< "This may lead to non-repeatable behaviour" << endl;
srand48(seed_);
osRandomSeed(seed_);
}
else if (reset)
{

73
src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.C
View File

@ -59,10 +59,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(p, iF),
phiName_("phi"),
rhoName_("rho"),
p0_(p.size(), 0.0),
totalPressureFvPatchScalarField(p, iF),
fanCurve_(),
direction_(ffdOut)
{}
@ -76,10 +73,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchScalarField(ptf, p, iF, mapper),
phiName_(ptf.phiName_),
rhoName_(ptf.rhoName_),
p0_(ptf.p0_, mapper),
totalPressureFvPatchScalarField(ptf, p, iF, mapper),
fanCurve_(ptf.fanCurve_),
direction_(ptf.direction_)
{}
@ -92,10 +86,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
const dictionary& dict
)
:
fixedValueFvPatchScalarField(p, iF),
phiName_(dict.lookupOrDefault<word>("phi", "phi")),
rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
p0_("p0", dict, p.size()),
totalPressureFvPatchScalarField(p, iF),
fanCurve_(dict),
direction_(fanFlowDirectionNames_.read(dict.lookup("direction")))
{
@ -109,10 +100,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
const fanPressureFvPatchScalarField& pfopsf
)
:
fixedValueFvPatchScalarField(pfopsf),
phiName_(pfopsf.phiName_),
rhoName_(pfopsf.rhoName_),
p0_(pfopsf.p0_),
totalPressureFvPatchScalarField(pfopsf),
fanCurve_(pfopsf.fanCurve_),
direction_(pfopsf.direction_)
{}
@ -124,10 +112,7 @@ Foam::fanPressureFvPatchScalarField::fanPressureFvPatchScalarField
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(pfopsf, iF),
phiName_(pfopsf.phiName_),
rhoName_(pfopsf.rhoName_),
p0_(pfopsf.p0_),
totalPressureFvPatchScalarField(pfopsf, iF),
fanCurve_(pfopsf.fanCurve_),
direction_(pfopsf.direction_)
{}
@ -144,7 +129,7 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs()
// Retrieve flux field
const surfaceScalarField& phi =
db().lookupObject<surfaceScalarField>(phiName_);
db().lookupObject<surfaceScalarField>(phiName());
const fvsPatchField<scalar>& phip =
patch().patchField<surfaceScalarField, scalar>(phi);
@ -161,7 +146,7 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs()
else if (phi.dimensions() == dimVelocity*dimArea*dimDensity)
{
const scalarField& rhop =
patch().lookupPatchField<volScalarField, scalar>(rhoName_);
patch().lookupPatchField<volScalarField, scalar>(rhoName());
aveFlowRate = dir*gSum(phip/rhop)/gSum(patch().magSf());
}
else
@ -174,51 +159,23 @@ void Foam::fanPressureFvPatchScalarField::updateCoeffs()
<< exit(FatalError);
}
// Normal flow through fan
if (aveFlowRate >= 0.0)
{
// Pressure drop for this flow rate
const scalar pdFan = fanCurve_(aveFlowRate);
// Pressure drop for this flow rate
const scalar pdFan = fanCurve_(max(aveFlowRate, 0.0));
operator==(p0_ - dir*pdFan);
}
// Reverse flow
else
{
// Assume that fan has stalled if flow reversed
// i.e. apply dp for zero flow rate
const scalar pdFan = fanCurve_(0);
// Flow speed across patch
scalarField Up = phip/(patch().magSf());
// Pressure drop associated withback flow = dynamic pressure
scalarField pdBackFlow = 0.5*magSqr(Up);
if (phi.dimensions() == dimVelocity*dimArea*dimDensity)
{
const scalarField& rhop =
patch().lookupPatchField<volScalarField, scalar>(rhoName_);
pdBackFlow /= rhop;
}
operator==(p0_ - dir*(pdBackFlow + pdFan));
}
fixedValueFvPatchScalarField::updateCoeffs();
totalPressureFvPatchScalarField::updateCoeffs
(
p0() - dir*pdFan,
patch().lookupPatchField<volVectorField, vector>(UName())
);
}
void Foam::fanPressureFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
os.writeKeyword("phi") << phiName_ << token::END_STATEMENT << nl;
os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
totalPressureFvPatchScalarField::write(os);
fanCurve_.write(os);
os.writeKeyword("direction")
<< fanFlowDirectionNames_[direction_] << token::END_STATEMENT << nl;
p0_.writeEntry("p0", os);
writeEntry("value", os);
}

20
src/finiteVolume/fields/fvPatchFields/derived/fanPressure/fanPressureFvPatchScalarField.H
View File

@ -25,7 +25,7 @@ Class
Foam::fanPressureFvPatchScalarField
Description
Assigns pressure inlet or outlet condition for a fan.
Assigns pressure inlet or outlet total pressure condition for a fan.
User specifies:
- pressure drop vs volumetric flow rate table (fan curve) file name;
@ -56,8 +56,8 @@ Description
\endverbatim
See Also
Foam::interpolationTable and
Foam::timeVaryingFlowRateInletVelocityFvPatchVectorField
Foam::totalPressureFvPatchScalarField and
Foam::interpolationTable
SourceFiles
fanPressureFvPatchScalarField.C
@ -67,8 +67,7 @@ SourceFiles
#ifndef fanPressureFvPatchScalarField_H
#define fanPressureFvPatchScalarField_H
#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"
#include "totalPressureFvPatchScalarField.H"
#include "interpolationTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -82,19 +81,10 @@ namespace Foam
class fanPressureFvPatchScalarField
:
public fixedValueFvPatchScalarField
public totalPressureFvPatchScalarField
{
// Private data
//- Name of the flux transporting the field
word phiName_;
//- Name of the density field
word rhoName_;
//- Total pressure
scalarField p0_;
//- Tabulated fan curve
interpolationTable<scalar> fanCurve_;

4
src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -112,7 +112,7 @@ void Foam::rotatingTotalPressureFvPatchScalarField::updateCoeffs()
+ rotationVelocity
);
totalPressureFvPatchScalarField::updateCoeffs(Up);
totalPressureFvPatchScalarField::updateCoeffs(p0(), Up);
}

2
src/finiteVolume/fields/fvPatchFields/derived/rotatingTotalPressure/rotatingTotalPressureFvPatchScalarField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License

22
src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -153,7 +153,11 @@ void Foam::totalPressureFvPatchScalarField::rmap
}
void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
void Foam::totalPressureFvPatchScalarField::updateCoeffs
(
const scalarField& p0p,
const vectorField& Up
)
{
if (updated())
{
@ -165,7 +169,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
if (psiName_ == "none" && rhoName_ == "none")
{
operator==(p0_ - 0.5*(1.0 - pos(phip))*magSqr(Up));
operator==(p0p - 0.5*(1.0 - pos(phip))*magSqr(Up));
}
else if (rhoName_ == "none")
{
@ -178,7 +182,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
operator==
(
p0_
p0p
/pow
(
(1.0 + 0.5*psip*gM1ByG*(1.0 - pos(phip))*magSqr(Up)),
@ -188,7 +192,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
}
else
{
operator==(p0_/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
operator==(p0p/(1.0 + 0.5*psip*(1.0 - pos(phip))*magSqr(Up)));
}
}
else if (psiName_ == "none")
@ -196,7 +200,7 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
const fvPatchField<scalar>& rho =
patch().lookupPatchField<volScalarField, scalar>(rhoName_);
operator==(p0_ - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
operator==(p0p - 0.5*rho*(1.0 - pos(phip))*magSqr(Up));
}
else
{
@ -220,7 +224,11 @@ void Foam::totalPressureFvPatchScalarField::updateCoeffs(const vectorField& Up)
void Foam::totalPressureFvPatchScalarField::updateCoeffs()
{
updateCoeffs(patch().lookupPatchField<volVectorField, vector>(UName_));
updateCoeffs
(
p0(),
patch().lookupPatchField<volVectorField, vector>(UName())
);
}

49
src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -157,6 +157,45 @@ public:
return UName_;
}
//- Return the name of the flux field
const word& phiName() const
{
return phiName_;
}
//- Return reference to the name of the flux field
// to allow adjustment
word& phiName()
{
return phiName_;
}
//- Return the name of the density field
const word& rhoName() const
{
return rhoName_;
}
//- Return reference to the name of the density field
// to allow adjustment
word& rhoName()
{
return rhoName_;
}
//- Return the name of the compressibility field
const word& psiName() const
{
return psiName_;
}
//- Return reference to the name of the compressibility field
// to allow adjustment
word& psiName()
{
return psiName_;
}
//- Return the heat capacity ratio
scalar gamma() const
{
@ -201,8 +240,12 @@ public:
// Evaluation functions
//- Update the coefficients associated with the patch field
// using the given patch velocity field
virtual void updateCoeffs(const vectorField& Up);
// using the given patch total pressure and velocity fields
virtual void updateCoeffs
(
const scalarField& p0p,
const vectorField& Up
);
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();

3
src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C
View File

@ -214,7 +214,8 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
rHat_AB
*(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
// Cohesion force
// Cohesion force, energy density multiplied by the area of
// particle-particle overlap
if (cohesion_)
{
fN_AB +=

42
src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -76,7 +76,8 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
typename CloudType::parcelType& p,
const point& site,
const WallSiteData<vector>& data,
scalar pREff
scalar pREff,
bool cohesion
) const
{
// wall patch index
@ -88,14 +89,18 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
scalar alpha = alpha_[wPI];
scalar b = b_[wPI];
scalar mu = mu_[wPI];
scalar cohesionEnergyDensity = cohesionEnergyDensity_[wPI];
cohesion = cohesion && cohesion_[wPI];
vector r_PW = p.position() - site;
vector U_PW = p.U() - data.wallData();
scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
scalar r_PW_mag = mag(r_PW);
vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;
@ -105,6 +110,16 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
rHat_PW
*(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));
// Cohesion force, energy density multiplied by the area of wall/particle
// overlap
if (cohesion)
{
fN_PW +=
-cohesionEnergyDensity
*mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
*rHat_PW;
}
p.f() += fN_PW;
vector USlip_PW =
@ -168,6 +183,8 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
alpha_(),
b_(),
mu_(),
cohesionEnergyDensity_(),
cohesion_(),
patchMap_(),
maxEstarIndex_(-1),
collisionResolutionSteps_
@ -212,6 +229,8 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
alpha_.setSize(nWallPatches);
b_.setSize(nWallPatches);
mu_.setSize(nWallPatches);
cohesionEnergyDensity_.setSize(nWallPatches);
cohesion_.setSize(nWallPatches);
scalar maxEstar = -GREAT;
@ -238,6 +257,13 @@ Foam::WallLocalSpringSliderDashpot<CloudType>::WallLocalSpringSliderDashpot
mu_[wPI] = readScalar(patchCoeffDict.lookup("mu"));
cohesionEnergyDensity_[wPI] = readScalar
(
patchCoeffDict.lookup("cohesionEnergyDensity")
);
cohesion_[wPI] = (mag(cohesionEnergyDensity_[wPI]) > VSMALL);
if (Estar_[wPI] > maxEstar)
{
maxEstarIndex_ = wPI;
@ -325,20 +351,22 @@ void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
p,
flatSitePoints[siteI],
flatSiteData[siteI],
pREff
pREff,
true
);
}
forAll(sharpSitePoints, siteI)
{
// Treating sharp sites like flat sites
// Treating sharp sites like flat sites, except suppress cohesion
evaluateWall
(
p,
sharpSitePoints[siteI],
sharpSiteData[siteI],
pREff
pREff,
false
);
}
}

9
src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallLocalSpringSliderDashpot/WallLocalSpringSliderDashpot.H
View File

@ -65,6 +65,12 @@ class WallLocalSpringSliderDashpot
//- Coefficient of friction in for tangential sliding
scalarList mu_;
//- Cohesion energy density [J/m^3]
scalarList cohesionEnergyDensity_;
// Switch cohesion on and off
boolList cohesion_;
//- Mapping the patch index to the model data
labelList patchMap_;
@ -115,7 +121,8 @@ class WallLocalSpringSliderDashpot
typename CloudType::parcelType& p,
const point& site,
const WallSiteData<vector>& data,
scalar pREff
scalar pREff,
bool cohesion
) const;

36
src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -77,16 +77,19 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
const point& site,
const WallSiteData<vector>& data,
scalar pREff,
scalar kN
scalar kN,
bool cohesion
) const
{
vector r_PW = p.position() - site;
vector U_PW = p.U() - data.wallData();
scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);
scalar r_PW_mag = mag(r_PW);
vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);
scalar normalOverlapMag = max(pREff - r_PW_mag, 0.0);
vector rHat_PW = r_PW/(r_PW_mag + VSMALL);
scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
@ -94,6 +97,16 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
rHat_PW
*(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
// Cohesion force, energy density multiplied by the area of wall/particle
// overlap
if (cohesion)
{
fN_PW +=
-cohesionEnergyDensity_
*mathematical::pi*(sqr(pREff) - sqr(r_PW_mag))
*rHat_PW;
}
p.f() += fN_PW;
vector USlip_PW =
@ -157,6 +170,11 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
alpha_(readScalar(this->coeffDict().lookup("alpha"))),
b_(readScalar(this->coeffDict().lookup("b"))),
mu_(readScalar(this->coeffDict().lookup("mu"))),
cohesionEnergyDensity_
(
readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
),
cohesion_(false),
collisionResolutionSteps_
(
readScalar
@ -183,6 +201,8 @@ Foam::WallSpringSliderDashpot<CloudType>::WallSpringSliderDashpot
Estar_ = 1/((1 - sqr(pNu))/pE + (1 - sqr(nu))/E);
Gstar_ = 1/(2*((2 + pNu - sqr(pNu))/pE + (2 + nu - sqr(nu))/E));
cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
}
@ -266,13 +286,14 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
flatSitePoints[siteI],
flatSiteData[siteI],
pREff,
kN
kN,
cohesion_
);
}
forAll(sharpSitePoints, siteI)
{
// Treating sharp sites like flat sites
// Treating sharp sites like flat sites, except suppress cohesion
evaluateWall
(
@ -280,7 +301,8 @@ void Foam::WallSpringSliderDashpot<CloudType>::evaluateWall
sharpSitePoints[siteI],
sharpSiteData[siteI],
pREff,
kN
kN,
false
);
}
}

9
src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/WallModel/WallSpringSliderDashpot/WallSpringSliderDashpot.H
View File

@ -65,6 +65,12 @@ class WallSpringSliderDashpot
//- Coefficient of friction in for tangential sliding
scalar mu_;
//- Cohesion energy density [J/m^3]
scalar cohesionEnergyDensity_;
// Switch cohesion on and off
bool cohesion_;
//- The number of steps over which to resolve the minimum
// harmonic approximation of the collision period
scalar collisionResolutionSteps_;
@ -110,7 +116,8 @@ class WallSpringSliderDashpot
const point& site,
const WallSiteData<vector>& data,
scalar pREff,
scalar kN
scalar kN,
bool cohesion
) const;

14
src/lagrangian/molecularDynamics/molecule/Make/files
View File

@ -1,8 +1,14 @@
constPropSite/constPropSite.C
clouds/baseClasses/moleculeCloud/moleculeCloud.C
polyatomic/polyatomic/polyatomic.C
polyatomic/polyatomic/polyatomicIO.C
molecules/constPropSite/constPropSite.C
polyatomic/polyatomicCloud/polyatomicCloud.C
molecules/monoatomic/monoatomic.C
molecules/monoatomic/monoatomicIO.C
molecules/polyatomic/polyatomic.C
molecules/polyatomic/polyatomicIO.C
/* controllers/basic/controllers/controllers.C
controllers/basic/stateController/stateController.C
controllers/basic/fluxController/fluxController.C */
LIB = $(FOAM_LIBBIN)/libmolecule

6
src/lagrangian/molecularDynamics/molecule/Make/options
View File

@ -2,13 +2,11 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
-I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude
LIB_LIBS = \
-lfiniteVolume \
-lmeshTools \
-llagrangian \
-lpotential \
-lmolecularMeasurements
-lpotential

372
src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.C → src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.C
View File

@ -23,36 +23,29 @@ License
\*----------------------------------------------------------------------------*/
#include "polyatomicCloud.H"
#include "MoleculeCloud.H"
#include "fvMesh.H"
#include "mathematicalConstants.H"
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTemplateTypeNameAndDebug(Cloud<polyatomic>, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::polyatomicCloud::buildConstProps()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::buildConstProps()
{
Info<< nl << "Reading polyatomicProperties dictionary." << endl;
const List<word>& idList(pot_.idList());
constPropList_.setSize(idList.size());
const List<word>& siteIdList(pot_.siteIdList());
IOdictionary polyatomicPropertiesDict
IOdictionary molPropertiesDict
(
IOobject
(
"polyatomicProperties",
this->name() + "Properties",
mesh_.time().constant(),
mesh_,
IOobject::MUST_READ_IF_MODIFIED,
@ -61,11 +54,13 @@ void Foam::polyatomicCloud::buildConstProps()
)
);
Info<< nl << "Reading " << molPropertiesDict.name() << endl;
forAll(idList, i)
{
const word& id = idList[i];
const dictionary& molDict(polyatomicPropertiesDict.subDict(id));
const dictionary& molDict(molPropertiesDict.subDict(id));
List<word> siteIdNames = molDict.lookup("siteIds");
@ -81,17 +76,14 @@ void Foam::polyatomicCloud::buildConstProps()
{
FatalErrorIn
(
"Foam::polyatomic::constantProperties::constantProperties"
"("
"const dictionary& dict"
")"
"void Foam::MoleculeCloud<MoleculeType>::buildConstProps()"
)
<< siteId << " site not found."
<< nl << abort(FatalError);
}
}
constPropList_[i] = polyatomic::constantProperties
constPropList_[i] = typename MoleculeType::constantProperties
(
molDict,
siteIds
@ -100,11 +92,15 @@ void Foam::polyatomicCloud::buildConstProps()
}
void Foam::polyatomicCloud::setSiteSizesAndPositions()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::setSiteSizesAndPositions()
{
forAllIter(polyatomicCloud, *this, mol)
forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
const polyatomic::constantProperties& cP = constProps(mol().id());
const typename MoleculeType::constantProperties& cP
(
constProps(mol().id())
);
mol().setSiteSizes(cP.nSites());
@ -113,14 +109,15 @@ void Foam::polyatomicCloud::setSiteSizesAndPositions()
}
void Foam::polyatomicCloud::buildCellOccupancy()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::buildCellOccupancy()
{
forAll(cellOccupancy_, cO)
{
cellOccupancy_[cO].clear();
}
forAllIter(polyatomicCloud, *this, mol)
forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
cellOccupancy_[mol().cell()].append(&mol());
}
@ -132,15 +129,16 @@ void Foam::polyatomicCloud::buildCellOccupancy()
}
void Foam::polyatomicCloud::calculatePairForce()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::calculatePairForce()
{
PstreamBuffers pBufs(Pstream::nonBlocking);
// Start sending referred data
il_.sendReferredData(cellOccupancy(), pBufs);
polyatomic* molI = NULL;
polyatomic* molJ = NULL;
MoleculeType* molI = NULL;
MoleculeType* molJ = NULL;
{
// Real-Real interactions
@ -155,7 +153,7 @@ void Foam::polyatomicCloud::calculatePairForce()
forAll(dil[d], interactingCells)
{
List<polyatomic*> cellJ =
List<MoleculeType*> cellJ =
cellOccupancy_[dil[d][interactingCells]];
forAll(cellJ, cellJMols)
@ -187,24 +185,24 @@ void Foam::polyatomicCloud::calculatePairForce()
const labelListList& ril = il_.ril();
List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
forAll(ril, r)
{
const List<label>& realCells = ril[r];
IDLList<polyatomic>& refMols = referredMols[r];
IDLList<MoleculeType>& refMols = referredMols[r];
forAllIter
(
IDLList<polyatomic>,
typename IDLList<MoleculeType>,
refMols,
refMol
)
{
forAll(realCells, rC)
{
List<polyatomic*> cellI = cellOccupancy_[realCells[rC]];
List<MoleculeType*> cellI = cellOccupancy_[realCells[rC]];
forAll(cellI, cellIMols)
{
@ -219,11 +217,12 @@ void Foam::polyatomicCloud::calculatePairForce()
}
void Foam::polyatomicCloud::calculateTetherForce()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::calculateTetherForce()
{
const tetherPotentialList& tetherPot(pot_.tetherPotentials());
forAllIter(polyatomicCloud, *this, mol)
forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
if (mol().tethered())
{
@ -246,16 +245,18 @@ void Foam::polyatomicCloud::calculateTetherForce()
}
void Foam::polyatomicCloud::calculateExternalForce()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::calculateExternalForce()
{
forAllIter(polyatomicCloud, *this, mol)
forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
mol().a() += pot_.gravity();
}
}
void Foam::polyatomicCloud::removeHighEnergyOverlaps()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::removeHighEnergyOverlaps()
{
Info<< nl << "Removing high energy overlaps, limit = "
<< pot_.potentialEnergyLimit()
@ -274,23 +275,23 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
// Real-Real interaction
polyatomic* molI = NULL;
polyatomic* molJ = NULL;
MoleculeType* molI = NULL;
MoleculeType* molJ = NULL;
{
DynamicList<polyatomic*> molsToDelete;
DynamicList<MoleculeType*> molsToDelete;
const labelListList& dil(il_.dil());
forAll(dil, d)
{
forAll(cellOccupancy_[d],cellIMols)
forAll(cellOccupancy_[d], cellIMols)
{
molI = cellOccupancy_[d][cellIMols];
forAll(dil[d], interactingCells)
{
List<polyatomic*> cellJ =
List<MoleculeType*> cellJ =
cellOccupancy_[dil[d][interactingCells]];
forAll(cellJ, cellJMols)
@ -322,35 +323,35 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
}
}
}
}
forAll(cellOccupancy_[d], cellIOtherMols)
{
molJ = cellOccupancy_[d][cellIOtherMols];
if (molJ > molI)
forAll(cellOccupancy_[d], cellIOtherMols)
{
if (evaluatePotentialLimit(*molI, *molJ))
molJ = cellOccupancy_[d][cellIOtherMols];
if (molJ > molI)
{
label idI = molI->id();
label idJ = molJ->id();
if
(
idI == idJ
|| findIndex(pot_.removalOrder(), idJ)
< findIndex(pot_.removalOrder(), idI)
)
if (evaluatePotentialLimit(*molI, *molJ))
{
if (findIndex(molsToDelete, molJ) == -1)
label idI = molI->id();
label idJ = molJ->id();
if
(
idI == idJ
|| findIndex(pot_.removalOrder(), idJ)
< findIndex(pot_.removalOrder(), idI)
)
{
molsToDelete.append(molJ);
if (findIndex(molsToDelete, molJ) == -1)
{
molsToDelete.append(molJ);
}
}
else if (findIndex(molsToDelete, molI) == -1)
{
molsToDelete.append(molI);
}
}
else if (findIndex(molsToDelete, molI) == -1)
{
molsToDelete.append(molI);
}
}
}
@ -376,19 +377,19 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
// Real-Referred interaction
{
DynamicList<polyatomic*> molsToDelete;
DynamicList<MoleculeType*> molsToDelete;
const labelListList& ril(il_.ril());
List<IDLList<polyatomic> >& referredMols = il_.referredParticles();
List<IDLList<MoleculeType> >& referredMols = il_.referredParticles();
forAll(ril, r)
{
IDLList<polyatomic>& refMols = referredMols[r];
IDLList<MoleculeType>& refMols = referredMols[r];
forAllIter
(
IDLList<polyatomic>,
typename IDLList<MoleculeType>,
refMols,
refMol
)
@ -401,7 +402,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
{
label cellI = realCells[rC];
List<polyatomic*> cellIMols = cellOccupancy_[cellI];
List<MoleculeType*> cellIMols = cellOccupancy_[cellI];
forAll(cellIMols, cIM)
{
@ -430,7 +431,7 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
== findIndex(pot_.removalOrder(), idJ)
)
{
// Remove one of the polyatomics
// Remove one of the molecules
// arbitrarily, assuring that a
// consistent decision is made for
// both real-referred pairs.
@ -470,17 +471,18 @@ void Foam::polyatomicCloud::removeHighEnergyOverlaps()
reduce(molsRemoved, sumOp<label>());
}
Info<< tab << molsRemoved << " polyatomics removed" << endl;
Info<< tab << molsRemoved << " molecules removed" << endl;
}
void Foam::polyatomicCloud::initialiseMolecules
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules
(
const IOdictionary& mdInitialiseDict
const dictionary& mdInitialiseDict
)
{
Info<< nl
<< "Initialising polyatomics in each zone specified in "
<< "Initialising molecules in each zone specified in "
<< mdInitialiseDict.name()
<< endl;
@ -488,7 +490,10 @@ void Foam::polyatomicCloud::initialiseMolecules
if (!cellZones.size())
{
FatalErrorIn("void Foam::polyatomicCloud::initialiseMolecules")
FatalErrorIn
(
"void Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
)
<< "No cellZones found in the mesh."
<< abort(FatalError);
}
@ -528,7 +533,10 @@ void Foam::polyatomicCloud::initialiseMolecules
if (latticeIds.size() != latticePositions.size())
{
FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
FatalErrorIn
(
"Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
)
<< "latticeIds and latticePositions must be the same "
<< " size." << nl
<< abort(FatalError);
@ -550,7 +558,10 @@ void Foam::polyatomicCloud::initialiseMolecules
if (numberDensity < VSMALL)
{
WarningIn("polyatomicCloud::initialiseMolecules")
WarningIn
(
"MoleculeCloud<MoleculeType>::initialiseMolecules"
)
<< "numberDensity too small, not filling zone "
<< zone.name() << endl;
@ -571,7 +582,7 @@ void Foam::polyatomicCloud::initialiseMolecules
{
label id = findIndex(pot_.idList(), latticeIds[i]);
const polyatomic::constantProperties& cP
const typename MoleculeType::constantProperties& cP
(
constProps(id)
);
@ -586,14 +597,16 @@ void Foam::polyatomicCloud::initialiseMolecules
if (massDensity < VSMALL)
{
WarningIn("polyatomicCloud::initialiseMolecules")
WarningIn
(
"MoleculeCloud<MoleculeType>::initialiseMolecules"
)
<< "massDensity too small, not filling zone "
<< zone.name() << endl;
continue;
}
latticeCellScale = pow
(
unitCellMass/(det(latticeCellShape)*massDensity),
@ -602,7 +615,10 @@ void Foam::polyatomicCloud::initialiseMolecules
}
else
{
FatalErrorIn("Foam::polyatomicCloud::initialiseMolecules")
FatalErrorIn
(
"Foam::MoleculeCloud<MoleculeType>::initialiseMolecules"
)
<< "massDensity or numberDensity not specified " << nl
<< abort(FatalError);
}
@ -697,8 +713,8 @@ void Foam::polyatomicCloud::initialiseMolecules
anchor += (R & (latticeCellShape & latticeAnchor));
// Continue trying to place polyatomics as long as at
// least one polyatomic is placed in each iteration.
// Continue trying to place molecule as long as at
// least one molecule is placed in each iteration.
// The "|| totalZoneMols == 0" condition means that the
// algorithm will continue if the origin is outside the
// zone.
@ -719,10 +735,6 @@ void Foam::polyatomicCloud::initialiseMolecules
bool partOfLayerInBounds = false;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// start of placement of polyatomics
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (n == 0)
{
// Special treatment is required for the first position,
@ -959,10 +971,6 @@ void Foam::polyatomicCloud::initialiseMolecules
}
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// end of placement of polyatomics
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if
(
totalZoneMols == 0
@ -971,11 +979,12 @@ void Foam::polyatomicCloud::initialiseMolecules
{
WarningIn
(
"Foam::polyatomicCloud::initialiseMolecules()"
"Foam::MoleculeCloud<MoleculeType>::"
"initialiseMolecules()"
)
<< "A whole layer of unit cells was placed "
<< "outside the bounds of the mesh, but no "
<< "polyatomics have been placed in zone '"
<< "molecules have been placed in zone '"
<< zone.name()
<< "'. This is likely to be because the zone "
<< "has few cells ("
@ -1001,7 +1010,8 @@ void Foam::polyatomicCloud::initialiseMolecules
}
void Foam::polyatomicCloud::createMolecule
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::createMolecule
(
const point& position,
label cell,
@ -1020,7 +1030,7 @@ void Foam::polyatomicCloud::createMolecule
if (cell == -1)
{
FatalErrorIn("Foam::polyatomicCloud::createMolecule")
FatalErrorIn("Foam::MoleculeCloud<MoleculeType>::createMolecule")
<< "Position specified does not correspond to a mesh cell." << nl
<< abort(FatalError);
}
@ -1033,59 +1043,31 @@ void Foam::polyatomicCloud::createMolecule
{
specialPosition = position;
special = polyatomic::SPECIAL_TETHERED;
special = MoleculeType::SPECIAL_TETHERED;
}
const polyatomic::constantProperties& cP(constProps(id));
const typename MoleculeType::constantProperties& cP(constProps(id));
vector v = equipartitionLinearVelocity(temperature, cP.mass());
v += bulkVelocity;
vector pi = vector::zero;
tensor Q = I;
if (!cP.pointMolecule())
{
pi = equipartitionAngularMomentum(temperature, cP);
scalar phi(rndGen_.scalar01()*twoPi);
scalar theta(rndGen_.scalar01()*twoPi);
scalar psi(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)
);
}
typename MoleculeType::trackingData td
(
*this,
MoleculeType::trackingData::tpAccess
);
addParticle
(
new polyatomic
new MoleculeType
(
mesh_,
position,
cell,
tetFace,
tetPt,
Q,
v,
vector::zero,
pi,
vector::zero,
temperature,
bulkVelocity,
specialPosition,
constProps(id),
cP,
td,
special,
id
)
@ -1093,11 +1075,12 @@ void Foam::polyatomicCloud::createMolecule
}
Foam::label Foam::polyatomicCloud::nSites() const
template<class MoleculeType>
Foam::label Foam::MoleculeCloud<MoleculeType>::nSites() const
{
label n = 0;
forAllConstIter(polyatomicCloud, *this, mol)
forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
n += constProps(mol().id()).nSites();
}
@ -1108,24 +1091,27 @@ Foam::label Foam::polyatomicCloud::nSites() const
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polyatomicCloud::polyatomicCloud
template<class MoleculeType>
Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
(
const word& cloudName,
const polyMesh& mesh,
const potential& pot,
bool readFields
)
:
Cloud<polyatomic>(mesh, "polyatomicCloud", false),
Cloud<MoleculeType>(mesh, cloudName, false),
moleculeCloud(),
mesh_(mesh),
pot_(pot),
cellOccupancy_(mesh_.nCells()),
il_(mesh_, pot_.pairPotentials().rCutMax(), false),
constPropList_(),
rndGen_(clock::getTime())
rndGen_(label(971501) + 1526*Pstream::myProcNo())
{
if (readFields)
{
polyatomic::readFields(*this);
MoleculeType::readFields(*this);
}
buildConstProps();
@ -1138,27 +1124,30 @@ Foam::polyatomicCloud::polyatomicCloud
}
Foam::polyatomicCloud::polyatomicCloud
template<class MoleculeType>
Foam::MoleculeCloud<MoleculeType>::MoleculeCloud
(
const word& cloudName,
const polyMesh& mesh,
const potential& pot,
const IOdictionary& mdInitialiseDict,
const dictionary& mdInitialiseDict,
bool readFields
)
:
Cloud<polyatomic>(mesh, "polyatomicCloud", false),
Cloud<MoleculeType>(mesh, cloudName, false),
moleculeCloud(),
mesh_(mesh),
pot_(pot),
il_(mesh_, 0.0, false),
il_(mesh_),
constPropList_(),
rndGen_(clock::getTime())
rndGen_(label(971501) + 1526*Pstream::myProcNo())
{
if (readFields)
{
polyatomic::readFields(*this);
MoleculeType::readFields(*this);
}
clear();
this->clear();
buildConstProps();
@ -1168,30 +1157,50 @@ Foam::polyatomicCloud::polyatomicCloud
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::polyatomicCloud::evolve()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::evolve()
{
polyatomic::trackingData td0(*this, 0);
Cloud<polyatomic>::move(td0, mesh_.time().deltaTValue());
typename MoleculeType::trackingData td0
(
*this,
MoleculeType::trackingData::tpFirstVelocityHalfStep
);
Cloud<MoleculeType>::move(td0, mesh_.time().deltaTValue());
polyatomic::trackingData td1(*this, 1);
Cloud<polyatomic>::move(td1, mesh_.time().deltaTValue());
typename MoleculeType::trackingData td1
(
*this,
MoleculeType::trackingData::tpLinearTrack
);
Cloud<MoleculeType>::move(td1, mesh_.time().deltaTValue());
polyatomic::trackingData td2(*this, 2);
Cloud<polyatomic>::move(td2, mesh_.time().deltaTValue());
typename MoleculeType::trackingData td2
(
*this,
MoleculeType::trackingData::tpRotationalTrack
);
Cloud<MoleculeType>::move(td2, mesh_.time().deltaTValue());
calculateForce();
polyatomic::trackingData td3(*this, 3);
Cloud<polyatomic>::move(td3, mesh_.time().deltaTValue());
typename MoleculeType::trackingData td3
(
*this,
MoleculeType::trackingData::tpSecondVelocityHalfStep
);
Cloud<MoleculeType>::move(td3, mesh_.time().deltaTValue());
info();
}
void Foam::polyatomicCloud::calculateForce()
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::calculateForce()
{
buildCellOccupancy();
// Set accumulated quantities to zero
forAllIter(polyatomicCloud, *this, mol)
forAllIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
mol().siteForces() = vector::zero;
@ -1208,45 +1217,36 @@ void Foam::polyatomicCloud::calculateForce()
}
void Foam::polyatomicCloud::applyConstraintsAndThermostats
(
const scalar targetTemperature,
const scalar measuredTemperature
)
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::info()
{
scalar temperatureCorrectionFactor =
sqrt(targetTemperature/max(VSMALL, measuredTemperature));
// Calculates and prints the mean momentum and energy in the system
// and the number of molecules.
Info<< "----------------------------------------" << nl
<< "Temperature equilibration" << nl
<< "Target temperature = "
<< targetTemperature << nl
<< "Measured temperature = "
<< measuredTemperature << nl
<< "Temperature correction factor = "
<< temperatureCorrectionFactor << nl
<< "----------------------------------------"
<< endl;
typename MoleculeType::trackingData td
(
*this,
MoleculeType::trackingData::tpAccess
);
forAllIter(polyatomicCloud, *this, mol)
{
mol().v() *= temperatureCorrectionFactor;
mol().pi() *= temperatureCorrectionFactor;
}
MoleculeType::info(td);
}
void Foam::polyatomicCloud::writeXYZ(const fileName& fName) const
template<class MoleculeType>
void Foam::MoleculeCloud<MoleculeType>::writeXYZ(const fileName& fName) const
{
OFstream os(fName);
os << nSites() << nl
<< "polyatomicCloud site positions in angstroms" << nl;
<< "MoleculeCloud<MoleculeType> site positions in angstroms" << nl;
forAllConstIter(polyatomicCloud, *this, mol)
forAllConstIter(typename MoleculeCloud<MoleculeType>, *this, mol)
{
const polyatomic::constantProperties& cP = constProps(mol().id());
const typename MoleculeType::constantProperties& cP
(
constProps(mol().id())
);
forAll(mol().sitePositions(), i)
{

134
src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloud.H → src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloud.H
View File

@ -22,22 +22,22 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::polyatomicCloud
Foam::MoleculeCloud
Description
SourceFiles
polyatomicCloudI.H
polyatomicCloud.C
MoleculeCloudI.H
MoleculeCloud.C
\*---------------------------------------------------------------------------*/
#ifndef polyatomicCloud_H
#define polyatomicCloud_H
#ifndef MoleculeCloud_H
#define MoleculeCloud_H
#include "Cloud.H"
#include "polyatomic.H"
#include "moleculeCloud.H"
#include "IOdictionary.H"
#include "potential.H"
#include "InteractionLists.H"
@ -51,65 +51,80 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class polyatomicCloud Declaration
Class MoleculeCloud Declaration
\*---------------------------------------------------------------------------*/
class polyatomicCloud
template<class MoleculeType>
class MoleculeCloud
:
public Cloud<polyatomic>
public Cloud<MoleculeType>,
public moleculeCloud
{
private:
// Private data
//-
const polyMesh& mesh_;
//-
const potential& pot_;
List<DynamicList<polyatomic*> > cellOccupancy_;
//-
List<DynamicList<MoleculeType*> > cellOccupancy_;
InteractionLists<polyatomic> il_;
//-
InteractionLists<MoleculeType> il_;
List<polyatomic::constantProperties> constPropList_;
//-
List<typename MoleculeType::constantProperties> constPropList_;
//-
Random rndGen_;
// Private Member Functions
//-
void buildConstProps();
//-
void setSiteSizesAndPositions();
//- Determine which polyatomics are in which cells
//- Determine which molecules are in which cells
void buildCellOccupancy();
//-
void calculatePairForce();
//-
inline void evaluatePair
(
polyatomic& molI,
polyatomic& molJ
MoleculeType& molI,
MoleculeType& molJ
);
//-
inline bool evaluatePotentialLimit
(
polyatomic& molI,
polyatomic& molJ
MoleculeType& molI,
MoleculeType& molJ
) const;
//-
void calculateTetherForce();
//-
void calculateExternalForce();
//-
void removeHighEnergyOverlaps();
void initialiseMolecules
(
const IOdictionary& mdInitialiseDict
);
//-
void initialiseMolecules(const dictionary& mdInitialiseDict);
//-
void createMolecule
(
const point& position,
@ -122,25 +137,14 @@ private:
const vector& bulkVelocity
);
//-
label nSites() const;
inline vector equipartitionLinearVelocity
(
scalar temperature,
scalar mass
);
inline vector equipartitionAngularMomentum
(
scalar temperature,
const polyatomic::constantProperties& cP
);
//- Disallow default bitwise copy construct
polyatomicCloud(const polyatomicCloud&);
MoleculeCloud(const MoleculeCloud&);
//- Disallow default bitwise assignment
void operator=(const polyatomicCloud&);
void operator=(const MoleculeCloud&);
public:
@ -148,59 +152,81 @@ public:
// Constructors
//- Construct given mesh and potential references
polyatomicCloud
MoleculeCloud
(
const word& cloudName,
const polyMesh& mesh,
const potential& pot,
bool readFields = true
);
//- Construct given mesh, potential and mdInitialiseDict
polyatomicCloud
MoleculeCloud
(
const word& cloudName,
const polyMesh& mesh,
const potential& pot,
const IOdictionary& mdInitialiseDict,
const dictionary& mdInitialiseDict,
bool readFields = true
);
// Member Functions
//- Evolve the polyatomics (move, calculate forces, control state etc)
//- Evolve the molecules (move, calculate forces, control state etc)
void evolve();
//-
void calculateForce();
void applyConstraintsAndThermostats
(
const scalar targetTemperature,
const scalar measuredTemperature
);
//- Print cloud information
void info();
// Access
//-
inline const polyMesh& mesh() const;
//-
inline const potential& pot() const;
inline const List<DynamicList<polyatomic*> >& cellOccupancy() const;
//-
inline const List<DynamicList<MoleculeType*> >&
cellOccupancy() const;
inline const InteractionLists<polyatomic>& il() const;
//-
inline const InteractionLists<MoleculeType>& il() const;
inline const List<polyatomic::constantProperties>
//-
inline const List<typename MoleculeType::constantProperties>
constProps() const;
inline const polyatomic::constantProperties&
constProps(label id) 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 polyatomic sites in XYZ format
//- Write molecule sites in XYZ format
void writeXYZ(const fileName& fName) const;
};
@ -211,7 +237,13 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "polyatomicCloudI.H"
#include "MoleculeCloudI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "MoleculeCloud.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

167
src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomicCloud/polyatomicCloudI.H → src/lagrangian/molecularDynamics/molecule/clouds/Templates/MoleculeCloud/MoleculeCloudI.H
View File

@ -29,10 +29,11 @@ using namespace Foam::constant;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline void Foam::polyatomicCloud::evaluatePair
template<class MoleculeType>
inline void Foam::MoleculeCloud<MoleculeType>::evaluatePair
(
polyatomic& molI,
polyatomic& molJ
MoleculeType& molI,
MoleculeType& molJ
)
{
const pairPotentialList& pairPot = pot_.pairPotentials();
@ -43,9 +44,15 @@ inline void Foam::polyatomicCloud::evaluatePair
label idJ = molJ.id();
const polyatomic::constantProperties& constPropI(constProps(idI));
const typename MoleculeType::constantProperties& constPropI
(
constProps(idI)
);
const polyatomic::constantProperties& constPropJ(constProps(idJ));
const typename MoleculeType::constantProperties& constPropJ
(
constProps(idJ)
);
forAll(constPropI.pairPotSites(), pI)
{
@ -149,10 +156,11 @@ inline void Foam::polyatomicCloud::evaluatePair
}
inline bool Foam::polyatomicCloud::evaluatePotentialLimit
template<class MoleculeType>
inline bool Foam::MoleculeCloud<MoleculeType>::evaluatePotentialLimit
(
polyatomic& molI,
polyatomic& molJ
MoleculeType& molI,
MoleculeType& molJ
) const
{
const pairPotentialList& pairPot = pot_.pairPotentials();
@ -163,9 +171,15 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
label idJ = molJ.id();
const polyatomic::constantProperties& constPropI(constProps(idI));
const typename MoleculeType::constantProperties& constPropI
(
constProps(idI)
);
const polyatomic::constantProperties& constPropJ(constProps(idJ));
const typename MoleculeType::constantProperties& constPropJ
(
constProps(idJ)
);
forAll(constPropI.pairPotSites(), pI)
{
@ -194,14 +208,18 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
if (rsIsJMag < SMALL)
{
WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
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 polyatomics"
<< " twice. Removing one of the polyatomics."
<< " parallel or a block filled with moleculess"
<< " twice. Removing one of the molecules."
<< endl;
return true;
@ -250,14 +268,18 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
if (rsIsJMag < SMALL)
{
WarningIn("polyatomicCloud::removeHighEnergyOverlaps()")
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 polyatomics"
<< " twice. Removing one of the polyatomics."
<< " parallel or a block filled with molecules"
<< " twice. Removing one of the molecules."
<< endl;
return true;
@ -288,7 +310,64 @@ inline bool Foam::polyatomicCloud::evaluatePotentialLimit
}
inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
// * * * * * * * * * * * * * * * 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
@ -303,10 +382,12 @@ inline Foam::vector Foam::polyatomicCloud::equipartitionLinearVelocity
}
inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
template<class MoleculeType>
inline Foam::vector
Foam::MoleculeCloud<MoleculeType>::equipartitionAngularMomentum
(
scalar temperature,
const polyatomic::constantProperties& cP
const typename MoleculeType::constantProperties& cP
)
{
scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
@ -332,52 +413,4 @@ inline Foam::vector Foam::polyatomicCloud::equipartitionAngularMomentum
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::polyMesh& Foam::polyatomicCloud::mesh() const
{
return mesh_;
}
inline const Foam::potential& Foam::polyatomicCloud::pot() const
{
return pot_;
}
inline const Foam::List<Foam::DynamicList<Foam::polyatomic*> >&
Foam::polyatomicCloud::cellOccupancy() const
{
return cellOccupancy_;
}
inline const Foam::InteractionLists<Foam::polyatomic>&
Foam::polyatomicCloud::il() const
{
return il_;
}
inline const Foam::List<Foam::polyatomic::constantProperties>
Foam::polyatomicCloud::constProps() const
{
return constPropList_;
}
inline const Foam::polyatomic::constantProperties&
Foam::polyatomicCloud::constProps(label id) const
{
return constPropList_[id];
}
inline Foam::Random& Foam::polyatomicCloud::rndGen()
{
return rndGen_;
}
// ************************************************************************* //

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

@ -0,0 +1,48 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,83 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,52 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,573 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,163 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,64 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,524 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,272 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,490 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,270 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
// ************************************************************************* //

52
src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.C → src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.C
View File

@ -21,12 +21,7 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
constPropSite
Description
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "constPropSite.H"
@ -71,6 +66,49 @@ 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;
}
// ************************************************************************* //

27
src/lagrangian/molecularDynamics/molecule/constPropSite/constPropSite.H → src/lagrangian/molecularDynamics/molecule/molecules/constPropSite/constPropSite.H
View File

@ -20,11 +20,13 @@ License
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
@ -35,14 +37,23 @@ SourceFiles
#define constPropSite_H
#include "vector.H"
#include "IFstream.H"
#include "OFstream.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
\*---------------------------------------------------------------------------*/
@ -144,6 +155,18 @@ public:
//-
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&);
};

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

@ -21,11 +21,6 @@ License
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
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -114,4 +109,31 @@ inline bool& Foam::constPropSite::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 Normal file
View File

@ -0,0 +1,199 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,418 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,328 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 Normal file
View File

@ -0,0 +1,305 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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;
}
// ************************************************************************* //

20
src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.C → src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.C
View File

@ -23,11 +23,19 @@ License
\*----------------------------------------------------------------------------*/
#include "polyatomicCloud.H"
#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
@ -81,7 +89,7 @@ bool Foam::polyatomic::move
const constantProperties& constProps(td.cloud().constProps(id_));
if (td.part() == 0)
if (td.part() == trackingData::tpFirstVelocityHalfStep)
{
// First leapfrog velocity adjust part, required before tracking+force
// part
@ -90,7 +98,7 @@ bool Foam::polyatomic::move
pi_ += 0.5*trackTime*tau_;
}
else if (td.part() == 1)
else if (td.part() == trackingData::tpLinearTrack)
{
// Leapfrog tracking part
@ -108,7 +116,7 @@ bool Foam::polyatomic::move
stepFraction() = 1.0 - tEnd/trackTime;
}
}
else if (td.part() == 2)
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
@ -149,7 +157,7 @@ bool Foam::polyatomic::move
setSitePositions(constProps);
}
else if (td.part() == 3)
else if (td.part() == trackingData::tpSecondVelocityHalfStep)
{
// Second leapfrog velocity adjust part, required after tracking+force
// part
@ -285,7 +293,7 @@ void Foam::polyatomic::hitWallPatch
)
{
// Use of the normal from tetIs is not required as
// hasWallImpactDistance for a polyatomicCloud is false.
// hasWallImpactDistance is false.
vector nw = normal();
nw /= mag(nw);

121
src/lagrangian/molecularDynamics/molecule/polyatomic/polyatomic/polyatomic.H → src/lagrangian/molecularDynamics/molecule/molecules/polyatomic/polyatomic.H
View File

@ -42,15 +42,13 @@ SourceFiles
#include "autoPtr.H"
#include "diagTensor.H"
#include "constPropSite.H"
#include "MoleculeCloud.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class polyatomicCloud;
/*---------------------------------------------------------------------------*\
Class polyatomic Declaration
\*---------------------------------------------------------------------------*/
@ -91,22 +89,19 @@ public:
//- Moment of intertia (in principal axis configiration)
diagTensor momentOfInertia_;
//-
scalar mass_;
// Private Member Functions
void setInteracionSiteBools
(
const List<word>& siteIds,
const List<word>& pairPotSiteIds
);
//-
bool linearMoleculeTest() const;
public:
//-
inline constantProperties();
//- Construct from dictionary
@ -118,22 +113,31 @@ public:
// 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;
};
@ -141,8 +145,22 @@ public:
//- Class used to pass tracking data to the trackToFace function
class trackingData
:
public particle::TrackingData<polyatomicCloud>
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_;
@ -151,9 +169,9 @@ public:
// Constructors
trackingData(polyatomicCloud& cloud, label part)
trackingData(MoleculeCloud<polyatomic>& cloud, trackPart part)
:
particle::TrackingData<polyatomicCloud>(cloud),
particle::TrackingData<MoleculeCloud<polyatomic> >(cloud),
part_(part)
{}
@ -176,10 +194,10 @@ private:
// bodyLocalVector = Q_.T() & globalVector
tensor Q_;
// Linear velocity of polyatomic
//- Linear velocity of polyatomic
vector v_;
// Total linear acceleration of polyatomic
//- Total linear acceleration of polyatomic
vector a_;
//- Angular momentum of polyatomic, in body local reference frame
@ -188,8 +206,10 @@ private:
//- Total torque on polyatomic, in body local reference frame
vector tau_;
//-
vector specialPosition_;
//-
scalar potentialEnergy_;
// - r_ij f_ij, stress dyad
@ -198,31 +218,59 @@ private:
// // - 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 from components
//- 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,
@ -279,53 +327,91 @@ public:
// 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;
@ -367,10 +453,15 @@ public:
// 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

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

@ -57,11 +57,12 @@ inline Foam::polyatomic::constantProperties::constantProperties
List<word> pairPotentialIds(dict.lookup("pairPotentialSiteIds"));
sites_.setSize(siteReferencePositions.size());
if
(
( siteIdNames.size() != sites_.size() )
|| ( siteReferencePositions.size() != sites_.size() )
|| ( siteCharges.size() != sites_.size() )
(siteIdNames.size() != sites_.size())
|| (siteCharges.size() != sites_.size())
)
{
FatalErrorIn
@ -72,12 +73,13 @@ inline Foam::polyatomic::constantProperties::constantProperties
")"
)
<< "Sizes of site id, charge and "
<< "referencePositions are not the same: " << sites_.size()
<< nl << abort(FatalError);
<< "referencePositions are not the same: " << nl
<< siteIdNames << nl
<< siteReferencePositions << nl
<< siteCharges << nl
<< abort(FatalError);
}
sites_.setSize(siteIdNames.size());
electrostaticSites_.setSize(sites_.size(), -1);
pairPotSites_.setSize(sites_.size(), -1);
@ -107,7 +109,7 @@ inline Foam::polyatomic::constantProperties::constantProperties
electrostaticSites_[electrostaticI++] = sI;
}
if (sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
if (!sites_[sI].pairPotentialSite() && !sites_[sI].electrostaticSite())
{
WarningIn
(
@ -303,6 +305,66 @@ inline Foam::polyatomic::constantProperties::constantProperties
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,
@ -315,10 +377,9 @@ inline Foam::polyatomic::polyatomic
const vector& pi,
const vector& tau,
const vector& specialPosition,
const constantProperties& constProps,
const constantProperties& cP,
const label special,
const label id
)
:
particle(mesh, position, cellI, tetFaceI, tetPtI),
@ -332,10 +393,10 @@ inline Foam::polyatomic::polyatomic
rf_(tensor::zero),
special_(special),
id_(id),
siteForces_(constProps.nSites(), vector::zero),
sitePositions_(constProps.nSites())
siteForces_(cP.nSites(), vector::zero),
sitePositions_(cP.nSites())
{
setSitePositions(constProps);
setSitePositions(cP);
}
@ -413,7 +474,7 @@ inline bool Foam::polyatomic::constantProperties::linearMolecule() const
}
inline bool Foam::polyatomic::constantProperties::pointMolecule() const
inline bool Foam::polyatomic::constantProperties::pointMolecule() const
{
return (momentOfInertia_.zz() < 0);
}

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

@ -25,7 +25,7 @@ License
#include "polyatomic.H"
#include "IOstreams.H"
#include "polyatomicCloud.H"
#include "Cloud.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -92,7 +92,11 @@ Foam::polyatomic::polyatomic
is.check
(
"Foam::polyatomic::polyatomic"
"(const Cloud<polyatomic>& cloud, Foam::Istream&), bool"
"("
"const polyMesh& mesh,"
"Istream& is,"
"bool readFields"
")"
);
}
@ -134,7 +138,8 @@ void Foam::polyatomic::readFields(Cloud<polyatomic>& mC)
mC.checkFieldIOobject(mC, id);
label i = 0;
forAllIter(polyatomicCloud, mC, iter)
forAllIter(typename Cloud<polyatomic>, mC, iter)
{
polyatomic& mol = iter();
@ -203,7 +208,7 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
);
label i = 0;
forAllConstIter(polyatomicCloud, mC, iter)
forAllConstIter(typename Cloud<polyatomic>, mC, iter)
{
const polyatomic& mol = iter();
@ -241,17 +246,114 @@ void Foam::polyatomic::writeFields(const Cloud<polyatomic>& mC)
orientation1.write();
orientation2.write();
orientation3.write();
}
Info<< "writeFields " << mC.name() << endl;
if (isA<polyatomicCloud>(mC))
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 polyatomicCloud& m = dynamic_cast<const polyatomicCloud&>(mC);
const label molId = mol().id();
scalar molMass(td.cloud().constProps(molId).mass());
totalMass += molMass;
//CentreOfMass += mol().position()*molMass;
}
m.writeXYZ
// if (nMols)
// {
// CentreOfMass /= totalMass;
// }
forAllConstIter(typename Cloud<polyatomic>, td.cloud(), mol)
{
const label molId = mol().id();
const polyatomic::constantProperties cP
(
m.mesh().time().timePath()/cloud::prefix/"polyatomicCloud.xmol"
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;
}
}

95
src/lagrangian/molecularDynamics/potential/potential/potential.C
View File

@ -93,7 +93,10 @@ void Foam::potential::setSiteIdList(const dictionary& moleculePropertiesDict)
}
void Foam::potential::potential::readPotentialDict()
void Foam::potential::potential::readPotentialDict
(
const dictionary& moleculePropertiesDict
)
{
Info<< nl << "Reading potential dictionary:" << endl;
@ -111,31 +114,28 @@ void Foam::potential::potential::readPotentialDict()
idList_ = List<word>(idListDict.lookup("idList"));
setSiteIdList
(
IOdictionary
(
IOobject
(
"moleculeProperties",
mesh_.time().constant(),
mesh_,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
setSiteIdList(moleculePropertiesDict);
List<word> pairPotentialSiteIdList
(
SubList<word>(siteIdList_, nPairPotIds_)
);
Info<< nl << "Unique site ids found: " << siteIdList_
<< nl << "Site Ids requiring a pair potential: "
<< pairPotentialSiteIdList
<< endl;
Info<< nl << "Unique site ids found:";
forAll(siteIdList_, i)
{
Info<< " " << siteIdList_[i];
}
Info << nl << "Site Ids requiring a pair potential:";
forAll(pairPotentialSiteIdList, i)
{
Info<< " " << pairPotentialSiteIdList[i];
}
Info<< nl;
List<word> tetherSiteIdList(0);
@ -229,37 +229,25 @@ void Foam::potential::potential::readPotentialDict()
if (potentialDict.found("external"))
{
Info<< nl << "Reading external forces:" << endl;
Info<< nl << "Reading external forces: ";
const dictionary& externalDict = potentialDict.subDict("external");
// gravity
externalDict.readIfPresent("gravity", gravity_);
}
Info<< nl << tab << "gravity = " << gravity_ << endl;
Info<< "gravity = " << gravity_ << nl << endl;
}
}
void Foam::potential::potential::readMdInitialiseDict
(
const IOdictionary& mdInitialiseDict,
IOdictionary& idListDict
const dictionary& mdInitialiseDict,
const dictionary& moleculePropertiesDict,
dictionary& idListDict
)
{
IOdictionary moleculePropertiesDict
(
IOobject
(
"moleculeProperties",
mesh_.time().constant(),
mesh_,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
);
DynamicList<word> idList;
DynamicList<word> tetherSiteIdList;
@ -286,11 +274,12 @@ void Foam::potential::potential::readMdInitialiseDict
(
"potential::readMdInitialiseDict"
"("
"const IOdictionary&, "
"IOdictionary&"
"const dictionary&, "
"dictionary&"
")"
) << "Molecule type " << id
<< " not found in moleculeProperties dictionary." << nl
<< " not found in " << moleculePropertiesDict.name()
<< " dictionary." << nl
<< abort(FatalError);
}
@ -341,8 +330,8 @@ void Foam::potential::potential::readMdInitialiseDict
(
"potential::readMdInitialiseDict"
"("
"const IOdictionary&, "
"IOdictionary&"
"const dictionary&, "
"dictionary&"
")"
) << "Tether id " << tetherSiteId
<< " not found as a site of any molecule in zone." << nl
@ -364,24 +353,34 @@ void Foam::potential::potential::readMdInitialiseDict
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::potential::potential(const polyMesh& mesh)
Foam::potential::potential
(
const polyMesh& mesh,
const dictionary& moleculePropertiesDict
)
:
mesh_(mesh)
{
readPotentialDict();
readPotentialDict(moleculePropertiesDict);
}
Foam::potential::potential
(
const polyMesh& mesh,
const IOdictionary& mdInitialiseDict,
IOdictionary& idListDict
const dictionary& mdInitialiseDict,
const dictionary& moleculePropertiesDict,
dictionary& idListDict
)
:
mesh_(mesh)
{
readMdInitialiseDict(mdInitialiseDict, idListDict);
readMdInitialiseDict
(
mdInitialiseDict,
moleculePropertiesDict,
idListDict
);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

40
src/lagrangian/molecularDynamics/potential/potential/potential.H
View File

@ -55,35 +55,49 @@ class potential
{
// Private data
//-
const polyMesh& mesh_;
//-
List<word> idList_;
//-
List<word> siteIdList_;
//-
label nPairPotIds_;
//-
scalar potentialEnergyLimit_;
//-
labelList removalOrder_;
//-
pairPotentialList pairPotentials_;
//-
tetherPotentialList tetherPotentials_;
//-
vector gravity_;
// Private Member Functions
//-
void setSiteIdList(const dictionary& moleculePropertiesDict);
void readPotentialDict();
//-
void readPotentialDict(const dictionary& moleculePropertiesDict);
//-
void readMdInitialiseDict
(
const IOdictionary& mdInitialiseDict,
IOdictionary& idListDict
const dictionary& mdInitialiseDict,
const dictionary& moleculePropertiesDict,
dictionary& idListDict
);
//- Disallow default bitwise copy construct
@ -98,14 +112,19 @@ public:
// Constructors
//- Construct from mesh reference
potential(const polyMesh& mesh);
potential
(
const polyMesh& mesh,
const dictionary& moleculePropertiesDict
);
//- Construct from mdInitialiseDict
potential
(
const polyMesh& mesh,
const IOdictionary& mdInitialiseDict,
IOdictionary& idListDict
const dictionary& mdInitialiseDict,
const dictionary& moleculePropertiesDict,
dictionary& idListDict
);
@ -117,22 +136,31 @@ public:
// Access
//-
inline label nIds() const;
//-
inline const List<word>& idList() const;
//-
inline const List<word>& siteIdList() const;
//-
inline scalar potentialEnergyLimit() const;
//-
inline label nPairPotentials() const;
//-
inline const labelList& removalOrder() const;
//-
inline const pairPotentialList& pairPotentials() const;
//-
inline const tetherPotentialList& tetherPotentials() const;
//-
inline const vector& gravity() const;
};

16
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -160,9 +160,19 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
if (closestSpokeSurface == -1)
{
FatalErrorIn("conformalVoronoiMesh::requiredAlignment")
<< "No secondary surface hit found in spoke search."
<< nl << exit(FatalError);
WarningIn
(
"conformalVoronoiMesh::requiredAlignment"
"("
"const Foam::point& pt"
") const"
)
<< "No secondary surface hit found in spoke search "
<< "using " << s
<< " spokes, try increasing alignmentSearchSpokes."
<< endl;
return I;
}
// Auxiliary alignment generated by spoke intersection normal.

51
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
View File

@ -1110,7 +1110,7 @@ Foam::conformalVoronoiMesh::collapseFace
// Output face and collapse axis for visualisation
Info<< "# Aspect ratio = " << aspectRatio << nl
Pout<< "# Aspect ratio = " << aspectRatio << nl
<< "# inertia = " << J << nl
<< "# determinant = " << detJ << nl
<< "# eigenvalues = " << eigenValues(J) << nl
@ -1120,17 +1120,17 @@ Foam::conformalVoronoiMesh::collapseFace
forAll(f, fPtI)
{
meshTools::writeOBJ(Info, pts[f[fPtI]]);
meshTools::writeOBJ(Pout, pts[f[fPtI]]);
}
Info<< "f";
Pout<< "f";
forAll(f, fPtI)
{
Info << " " << fPtI + 1;
Pout << " " << fPtI + 1;
}
Info<< endl;
Pout<< endl;
return fcmNone;
}
@ -1402,7 +1402,7 @@ Foam::conformalVoronoiMesh::collapseFace
// {
// // Output face and collapse axis for visualisation
// Info<< "# Aspect ratio = " << aspectRatio << nl
// Pout<< "# Aspect ratio = " << aspectRatio << nl
// << "# determinant = " << detJ << nl
// << "# collapseAxis = " << collapseAxis << nl
// << "# mode = " << mode << nl
@ -1412,33 +1412,33 @@ Foam::conformalVoronoiMesh::collapseFace
// scalar scale = 2.0*mag(fC - pts[f[0]]);
// meshTools::writeOBJ(Info, fC);
// meshTools::writeOBJ(Info, fC + scale*collapseAxis);
// meshTools::writeOBJ(Pout, fC);
// meshTools::writeOBJ(Pout, fC + scale*collapseAxis);
// Info<< "f 1 2" << endl;
// Pout<< "f 1 2" << endl;
// forAll(f, fPtI)
// {
// meshTools::writeOBJ(Info, pts[f[fPtI]]);
// meshTools::writeOBJ(Pout, pts[f[fPtI]]);
// }
// Info<< "f";
// Pout<< "f";
// forAll(f, fPtI)
// {
// Info << " " << fPtI + 3;
// Pout << " " << fPtI + 3;
// }
// Info<< nl << "# " << d << endl;
// Pout<< nl << "# " << d << endl;
// Info<< "# " << d.first() << " " << d.last() << endl;
// Pout<< "# " << d.first() << " " << d.last() << endl;
// forAll(d, dI)
// {
// meshTools::writeOBJ(Info, fC + (d[dI] - dShift)*collapseAxis);
// meshTools::writeOBJ(Pout, fC + (d[dI] - dShift)*collapseAxis);
// }
// Info<< endl;
// Pout<< endl;
// }
return mode;
@ -1490,7 +1490,7 @@ void Foam::conformalVoronoiMesh::deferredCollapseFaceSet
}
}
Info<< "facesToCollapse" << nl << faceLabels << endl;
Pout<< "facesToCollapse" << nl << faceLabels << endl;
}
@ -1632,7 +1632,8 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
if (checkFaces.size() < fAreas.size())
{
Info<< "Excluding " << fAreas.size() - checkFaces.size()
Info<< "Excluding "
<< returnReduce(fAreas.size() - checkFaces.size(), sumOp<label>())
<< " faces from check, < " << faceAreaLimit << " area" << endl;
}
@ -1655,7 +1656,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
{
if (cells[cI].size() < 4 && cells[cI].size() > 0)
{
// Info<< "cell " << cI << " " << cells[cI]
// Pout<< "cell " << cI << " " << cells[cI]
// << " has " << cells[cI].size() << " faces."
// << endl;
@ -1669,7 +1670,8 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
}
Info<< " cells with more than 1 but fewer than 4 faces : "
<< nInvalidPolyhedra << endl;
<< returnReduce(nInvalidPolyhedra, sumOp<label>())
<< endl;
// Check for cells with one internal face only
@ -1712,7 +1714,8 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
}
Info<< " cells with with zero or one non-boundary face : "
<< oneInternalFaceCells << endl;
<< returnReduce(oneInternalFaceCells, sumOp<label>())
<< endl;
}
@ -1795,7 +1798,7 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::checkPolyMeshQuality
}
Info<< nl << "Maximum number of filter limits applied: "
<< maxFilterCount << endl;
<< returnReduce(maxFilterCount, maxOp<label>()) << endl;
return wrongFaces;
}
@ -2489,7 +2492,9 @@ void Foam::conformalVoronoiMesh::removeUnusedPoints
inplaceReorder(oldToNew, pts);
Info<< " Removing " << pts.size() - pointI << " unused points"
Info<< " Removing "
<< returnReduce(pts.size() - pointI, sumOp<label>())
<< " unused points"
<< endl;
pts.setSize(pointI);

5
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
View File

@ -696,9 +696,10 @@ void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
}
}
if (!protrudingCells.empty())
if (returnReduce(protrudingCells.size(), sumOp<label>()) > 0)
{
Pout<< nl << "Found " << protrudingCells.size()
Info<< nl << "Found "
<< returnReduce(protrudingCells.size(), sumOp<label>())
<< " cells protruding from the surface, writing cellSet "
<< protrudingCells.name()
<< endl;

17
src/postProcessing/functionObjects/IO/partialWrite/partialWrite.C
View File

@ -27,6 +27,7 @@ License
#include "dictionary.H"
#include "Time.H"
#include "IOobjectList.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -123,14 +124,26 @@ void Foam::partialWrite::write()
else
{
// Delete all but marked objects
fileName dbDir;
if (isA<polyMesh>(obr_))
{
dbDir = dynamic_cast<const polyMesh&>(obr_).dbDir();
}
IOobjectList objects(obr_, obr_.time().timeName());
forAllConstIter(HashPtrTable<IOobject>, objects, iter)
{
if (!objectNames_.found(iter()->name()))
{
const fileName f = obr_.time().timePath()/iter()->name();
//Pout<< " rm " << f << endl;
const fileName f =
obr_.time().timePath()
/dbDir
/iter()->name();
if (debug)
{
Pout<< " rm " << f << endl;
}
rm(f);
}
}

4
tutorials/combustion/fireFoam/les/oppositeBurningPanels/Allrun
View File

@ -32,4 +32,8 @@ decomposePar -region panelRegion > log.decomposParPanelRegion.log 2>&1
runParallel `getApplication` 6
paraFoam -touch
paraFoam -touch -region panelRegion
# ----------------------------------------------------------------- end-of-file

2
tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperEmptying/constant/kinematicCloudProperties
View File

@ -112,6 +112,7 @@ subModels
alpha 0.12;
b 1.5;
mu 0.43;
cohesionEnergyDensity 0;
}
frontAndBack
{
@ -120,6 +121,7 @@ subModels
alpha 0.12;
b 1.5;
mu 0.1;
cohesionEnergyDensity 0;
}
};
}

2
tutorials/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState/constant/kinematicCloudProperties
View File

@ -121,6 +121,7 @@ subModels
alpha 0.12;
b 1.5;
mu 0.43;
cohesionEnergyDensity 0;
}
frontAndBack
{
@ -129,6 +130,7 @@ subModels
alpha 0.12;
b 1.5;
mu 0.1;
cohesionEnergyDensity 0;
}
};
}

3
tutorials/lagrangian/reactingParcelFilmFoam/cylinder/Allrun.pre
View File

@ -8,3 +8,6 @@ runApplication setSet -batch wallFilmRegion.setSet
mv log.setSet log.wallFilmRegion.setSet
runApplication extrudeToRegionMesh -overwrite
paraFoam -touch
paraFoam -touch -region wallFilmRegion

3
tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/Allrun.pre
View File

@ -32,3 +32,6 @@ cp -r system/wallFilmRegion.org system/wallFilmRegion
find ./0 -maxdepth 1 -type f -exec \
sed -i "s/wallFilm/\"(region0_to.*)\"/g" {} \;
paraFoam -touch
paraFoam -touch -region wallFilmRegion

2
tutorials/lagrangian/reactingParcelFilmFoam/rivuletPanel/Allrun.pre
View File

@ -10,3 +10,5 @@ mv log.setSet log.wallFilmRegion.setSet
runApplication extrudeToRegionMesh -overwrite
paraFoam -touch
paraFoam -touch -region wallFilmRegion

3
tutorials/lagrangian/reactingParcelFilmFoam/splashPanel/Allrun.pre
View File

@ -14,3 +14,6 @@ runApplication setSet -region wallFilmRegion -batch createWallFilmRegionPatches.
mv log.setSet log.createWallFilmRegionPatches.setSet
runApplication createPatch -region wallFilmRegion -overwrite
paraFoam -touch
paraFoam -touch -region wallFilmRegion