Particle insertion

- particleInsertion is now transfered to time folders to track particle insertion progress
- collision check is added for new inserted particles and old particles in the simulaiton.
- minor changes to dataEntry and triple
This commit is contained in:
Hamidreza Norouzi 2024-04-15 13:19:03 -07:00
parent d6798b5097
commit 79b987c711
18 changed files with 397 additions and 167 deletions

View File

@ -34,7 +34,7 @@ pFlow::particleWallContactSearchs<method>::particleWallContactSearchs
domainBox_(domain),
updateInterval_
(
dict.getValOrSet<uint32>("updateInterval", 1)
max(dict.getValOrSet<uint32>("updateInterval", 1),1u)
)
{

View File

@ -18,31 +18,7 @@ Licence:
-----------------------------------------------------------------------------*/
/*template<typename ShapeType>
bool pFlow::Insertion<ShapeType>::readInsertionDict
(
const dictionary& dict
)
{
if(!insertion::readInsertionDict(dict)) return false;
regions_.clear();
if( !this->isActive() )
{
return true;
}
wordList regionDicNames = dict.dictionaryKeywords();
for(auto& name:regionDicNames)
{
REPORT(2)<<"reading insertion region "<< greenText(name)<<endREPORT;
regions_.push_backSafe(dict.subDict(name), shapes_);
}
return true;
}
template<typename ShapeType>
bool pFlow::Insertion<ShapeType>::writeInsertionDict
@ -50,7 +26,8 @@ bool pFlow::Insertion<ShapeType>::writeInsertionDict
dictionary& dict
)const
{
if( !insertion::writeInsertionDict(dict) ) return false;
if(!insertion::writeInsertionDict(dict))return false;
if( !this->isActive() ) return true;
@ -60,16 +37,18 @@ bool pFlow::Insertion<ShapeType>::writeInsertionDict
if( !regions_[i].write(rgnDict) )
{
fatalErrorInFunction<<
"Error in writing to dictionary "<<rgnDict.globalName()<<endl;
return false;
}
}
return true;
}*/
}
template<typename ShapeType>
bool
pFlow::Insertion<ShapeType>::setInsertionRegions()
pFlow::Insertion<ShapeType>::readInsertionDict()
{
regions_.clear();
@ -102,7 +81,7 @@ pFlow::Insertion<ShapeType>::Insertion(
insertion(prtcl),
shapes_(shapes)
{
if(!setInsertionRegions())
if(!readInsertionDict())
{
fatalErrorInFunction;
fatalExit;

View File

@ -62,12 +62,11 @@ private:
// - insertion regions
ListPtr<InsertionRegion<ShapeType>> regions_;
bool readInsertionDict();
bool setInsertionRegions();
protected:
// bool readInsertionDict(const dictionary& dict);
// bool writeInsertionDict(dictionary& dict)const;
bool writeInsertionDict(dictionary& dict)const override;
public:

View File

@ -55,8 +55,8 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
uint32 iter,
real t,
real dt,
wordVector& names,
realx3Vector& pos,
wordVector& names,
realx3Vector& positions,
bool& insertionOccured
)
{
@ -65,65 +65,185 @@ bool pFlow::InsertionRegion<ShapeType>::insertParticles
if(!insertionTime(iter, t, dt)) return true;
uint32 newNum = numberToBeInserted(iter, t, dt);
if(newNum == 0) return true;
// get the internal box
auto internalBox = pStruct().internalDomainBox();
if( !(internalBox.minPoint() < internalBox.maxPoint()))
{
WARNING<<"Minimum point of internal point is not lower than "<<
"the maximum point \n"<<
"minimum point: "<< internalBox.minPoint()<<
"\nmaximum point:"<<internalBox.maxPoint()<<END_WARNING;
return false;
}
names.reserve(max(newNum,names.capacity()));
pos.reserve(max(newNum,pos.capacity()));
names.clear();
pos.clear();
realVector diams("diams", newNum, 0, RESERVE());
size_t n = 0;
uint32 idx;
/// reserve enough space
names.reserve(max(newNum,names.capacity()));
names.clear();
positions.reserve(max(newNum,positions.capacity()));
positions.clear();
auto allPositions = realx3Vector("allPositions");
auto allDiameters = realVector("allDiameters");
auto& mix = mixture();
auto& pReg = pRegion();
real maxDiam = shapes_.maxBoundingSphere();
auto minP = pReg.minPoint() - maxDiam;
auto maxP = pReg.maxPoint() + maxDiam;
if(Insertion().checkForCollision())
{
// check for collision with already inserted particles
// so, use selector to select particles in the simulation
auto bDict = dictionary("boxInfo");
bDict.add("min", minP);
bDict.add("max", maxP);
auto selector = pStructSelector::create(
"box",
pStruct(),
bDict);
allPositions = selector().selectedPointPositions();
allDiameters = selectedFieldVals<real>
(
selector(),
Insertion().diameterName()
);
}
auto collCheck = collisionCheck(
{minP, maxP},
maxDiam,
allPositions,
allDiameters);
uint32 numInserted = 0;
uint32 idx;
word name = mix.getNextShapeName();
shapes_.shapeNameToIndex(name, idx);
real d = shapes_.boundingDiameter(idx);
for(uint32 i=0; i< 100*newNum ; ++i)
{
realx3 p = pReg.peek();
// check if point is inside internal box
if(!internalBox.isInside(p))continue;
if( !checkForContact(pos, diams, p, d) )
if( collCheck.checkPoint( p, d) )
{
names.push_back(name);
pos.push_back(p);
diams.push_back(d);
n++;
if( n == newNum )
{
addToNumInserted(newNum);
insertionOccured = true;
return true;
}
positions.push_back(p);
numInserted++;
if( numInserted == newNum ) break;
// add this new particle to collision check set
allDiameters.push_back(d);
allPositions.push_back(p);
collCheck.mapLastAddedParticle();
// obtain next shape name and diameter
name = mix.getNextShapeName();
shapes_.shapeNameToIndex(name, idx);
d = shapes_.boundingDiameter(idx);
}
}
addToNumInserted(n);
insertionOccured = true;
return false;
}
addToNumInserted(numInserted);
return numInserted == newNum;
}
/*if(!checkForCollision)
{
realVector diams("diams", newNum, 0, RESERVE());
uint32 idx;
word name = mix.getNextShapeName();
shapes_.shapeNameToIndex(name, idx);
real d = shapes_.boundingDiameter(idx);
for(uint32 i=0; i< 100*newNum ; ++i)
{
realx3 p = pReg.peek();
// check if point is inside internal box
if(!internalBox.isInside(p))continue;
if( !checkForContact(positions, diams, p, d) )
{
names.push_back(name);
positions.push_back(p);
diams.push_back(d);
numInserted++;
if( numInserted == newNum ) break;
name = mix.getNextShapeName();
shapes_.shapeNameToIndex(name, idx);
d = shapes_.boundingDiameter(idx);
}
}
}
else
{
real maxDiam = shapes_.maxBoundingSphere();
auto minP = pReg.minPoint() - maxDiam;
auto maxP = pReg.maxPoint() + maxDiam;
auto bDict = dictionary("boxInfo");
bDict.add("min", minP);
bDict.add("max", maxP);
auto selector = pStructSelector::create(
"box",
pStruct(),
bDict);
auto allPositions = selector().selectedPointPositions();
auto allDiameters = selectedFieldVals<real>(selector(), "diameter");
auto collCheck = collisionCheck(
{minP, maxP},
maxDiam,
allPositions,
allDiameters);
uint32 idx;
word name = mix.getNextShapeName();
shapes_.shapeNameToIndex(name, idx);
real d = shapes_.boundingDiameter(idx);
for(uint32 i=0; i< 100*newNum ; ++i)
{
realx3 p = pReg.peek();
// check if point is inside internal box
if(!internalBox.isInside(p))continue;
if( collCheck.checkPoint( p, d) )
{
names.push_back(name);
positions.push_back(p);
numInserted++;
if( numInserted == newNum ) break;
// add this new particle to collision check set
allDiameters.push_back(d);
allPositions.push_back(p);
collCheck.mapLastAddedParticle();
// obtain next shape name and diameter
name = mix.getNextShapeName();
shapes_.shapeNameToIndex(name, idx);
d = shapes_.boundingDiameter(idx);
}
}
}*/

View File

@ -21,9 +21,15 @@ Licence:
#ifndef __InsertionRegion_hpp__
#define __InsertionRegion_hpp__
#include "insertionRegion.hpp"
#include "dictionary.hpp"
#include "pointStructure.hpp"
#include "insertionRegion.hpp"
#include "insertion.hpp"
#include "collisionCheck.hpp"
#include "pStructSelector.hpp"
#include "fieldSelector.hpp"
namespace pFlow
{

View File

@ -49,16 +49,29 @@ pFlow::collisionCheck::checkPoint(const realx3& p, const real d) const
if(!searchBox_.isInside(p)) return false;
const auto ind = pointIndex(p);
uint32 n = head_(ind.x(), ind.y(), ind.z());
while( n != -1)
const auto startInd = max(ind - 1 ,int32x3(0));
const auto endInd = min( ind+1 ,nCells_-1);
for(int32 i=startInd.x(); i<=endInd.x(); i++)
{
if( (position_[n]-p).length() - 0.5*(diameters_[n]+d )<0.0 )
for(int32 j=startInd.y(); j<=endInd.y(); j++)
{
return false;
for(int32 k=startInd.z(); k<=endInd.z(); k++)
{
uint32 n = head_(i, j, k);
while( n != -1)
{
if( ((position_[n]-p).length() - 0.5*(diameters_[n]+d )) <= 0.0 )
{
return false;
}
n = next_[n];
}
}
}
n = next_[n];
}
return true;
}

View File

@ -30,13 +30,42 @@ pFlow::insertion::insertion(particles& prtcl)
insertionFile__,
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER
objectFile::WRITE_ALWAYS
),
&prtcl.control().caseSetup()
&prtcl.time()
),
particles_(prtcl)
{
readInsertionDict(*this);
// this means that insertion file exist in time folder
if( IOobject::implyRead() )
{
readFromFile_ = true;
} // look inside the caseSetup folder if it exist
else
{
// read dictionary from caseSetup folder
fileDictionary caseFile(
objectFile(
insertionFile__,
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_NEVER),
&prtcl.control().caseSetup());
// check if read happened
if(caseFile.implyRead())
{
readFromFile_ = true;
// assign it to this dictionary
fileDictionary::dictionary::operator=(caseFile);
}
}
if( readFromFile_)
{
readInsertionDict();
}
}
const pFlow::pointStructure&
@ -45,27 +74,18 @@ pFlow::insertion::pStruct() const
return particles_.pStruct();
}
bool
pFlow::insertion::read(iIstream& is, const IOPattern& iop)
{
if (fileDictionary::read(is, iop))
{
readFromFile_ = true;
return true;
}
else
{
return false;
}
}
bool
pFlow::insertion::readInsertionDict(const dictionary& dict)
pFlow::insertion::readInsertionDict()
{
active_ = dict.getVal<Logical>("active");
active_ = getVal<Logical>("active");
if (active_)
{
checkForCollision_ = getVal<Logical>("checkForCollision");
REPORT(1) << "Particle insertion mechanism is " << Yellow_Text("active")
<< " in the simulation." << END_REPORT;
}
@ -78,26 +98,56 @@ pFlow::insertion::readInsertionDict(const dictionary& dict)
return true;
}
/*
bool pFlow::insertion::writeInsertionDict
(
dictionary& dict
)const
bool
pFlow::insertion::writeInsertionDict(dictionary& dict) const
{
if(!dict.add("active", active_) )
{
fatalErrorInFunction<<
" error in writing active to dictionary
"<<dict.globalName()<<endl; return false;
}
if (!dict.add("active", active_))
{
fatalErrorInFunction <<"Error in writing active to dictionary "
<<dict.globalName()<<endl;
return false;
}
if(!dict.add("checkForCollision", checkForCollision_) )
{
fatalErrorInFunction<<
" error in writing checkForCollision to
dictionary
"<<dict.globalName()<<endl; return false;
}
if (!dict.add("checkForCollision", checkForCollision_))
{
fatalErrorInFunction <<
"Error in writing checkForCollision to dictionary "<<
dict.globalName()<<endl;
return false;
}
return true;
return true;
}
bool pFlow::insertion::write(iOstream & os, const IOPattern & iop) const
{
dictionary newDict(fileDictionary::dictionary::name(), true);
if( iop.thisProcWriteData() )
{
if( !writeInsertionDict(newDict) ||
!newDict.write(os))
{
fatalErrorInFunction<<
" error in writing to dictionary "<< newDict.globalName()<<endl;
return false;
}
}
return true;
}
/*bool
pFlow::insertion::read(iIstream& is, const IOPattern& iop)
{
if (fileDictionary::read(is, iop))
{
readFromFile_ = true;
return true;
}
else
{
return false;
}
}*/

View File

@ -46,16 +46,22 @@ private:
/// to be inserted due to collision
Logical increaseVelocity_ = "No";
word diameterName_ = "diameter";
word velocityName_ = "velocity";
/// Ref to particles
particles& particles_;
bool readFromFile_ = false;
/// Read from dictionary
bool readInsertionDict(const dictionary& dict);
bool readInsertionDict();
protected:
/// Write to dictionary
// bool writeInsertionDict(dictionary& dict)const;
virtual bool writeInsertionDict(dictionary& dict)const;
public:
@ -71,7 +77,7 @@ public:
/// is Insertion active
inline bool isActive() const
{
return active_();
return readFromFile_ && active_();
}
inline bool checkForCollision() const
@ -96,14 +102,24 @@ public:
return readFromFile_;
}
const word& diameterName()const
{
return diameterName_;
}
const word& velocityName()const
{
return velocityName_;
}
/// read from stream
bool read(iIstream& is, const IOPattern& iop) override;
//bool read(iIstream& is, const IOPattern& iop) override;
/*/// read from iIstream
virtual bool read(iIstream& is) = 0;
virtual bool read(iIstream& is) = 0;*/
/// write to iOstream
virtual bool write(iOstream& os)const = 0;*/
bool write(iOstream& os, const IOPattern& iop)const override ;
};
}

View File

@ -101,9 +101,15 @@ pFlow::insertionRegion::readInsertionRegion(const dictionary& dict)
bool
pFlow::insertionRegion::writeInsertionRegion(dictionary& dict) const
{
if (!dict.add("type", type_))
if(!dict.add("rate", rate_))
return false;
if(!tControl_.write(dict))
return false;
if (!dict.add("regionType", type_))
return false;
if (pRegion_)
{
auto& prDict = dict.subDictOrCreate(type_ + "Info");
@ -118,11 +124,11 @@ pFlow::insertionRegion::writeInsertionRegion(dictionary& dict) const
return false;
}
/*if(setFields_)
if(setFieldDict_)
{
auto& sfDict = dict.subDictOrCreate("setFields");
setFields_().write(sfDict);
}*/
if(!dict.addDict("setFields", setFieldDict_()))
return false;
}
return true;
}

View File

@ -159,6 +159,11 @@ public:
return type_;
}
const auto& dict()const
{
return dict_;
}
const auto& Insertion() const
{
return insertion_;
@ -206,16 +211,13 @@ public:
return false;
return readInsertionRegion(dict);
}
}*/
/// write to dictionary
bool write(dictionary& dict) const
{
if (!timeFlowControl::write(dict))
return false;
return writeInsertionRegion(dict);
}*/
}
};
} // pFlow

View File

@ -58,7 +58,9 @@ private:
uint32Vector numberInserted_{"numberInserted"};
/// Current number of inserted
uint32Vector current_{"currentInserted"};
uint32Vector current_{"currentInserted"};
uint32 lastPeaked_ = 0;
public:

View File

@ -61,7 +61,7 @@ pFlow::particleIdHandler::hearChanges(
auto idRange = getIdRange(numNew);
uint32Vector newId("newId",numNew,numNew,RESERVE());
fillSequence(newId, idRange.first);
output<< "id "<< idRange<<endl;
return this->field().insertSetElement(indices, newId);
}
else

View File

@ -26,8 +26,7 @@ bool pFlow::internalField<T, MemorySpace>::insert(const anyList& varList)
const auto& indices = varList.getObject<uint32IndexContainer>(
eventName);
bool success = false;
output<<"insert for field "<< name()<<endl;
if(varList.contains(name()))
{
// a single value is assigned

View File

@ -121,6 +121,9 @@ bool pFlow::dataEntry::writeDataEntry
}
else
{
if(tok != tokens.cbegin())os<<spaceToken();
os << *tok;
}
lastPuncBegin = true;
@ -141,6 +144,7 @@ bool pFlow::dataEntry::writeDataEntry
}
else
{
if(!lastPuncBegin&&applySpace) os<<spaceToken();
os << *tok;
lastPuncEnd = false;

View File

@ -27,6 +27,8 @@ pFlow::baseTimeControl::baseTimeControl
const word& intervalPrefix,
real defStartTime
)
:
intervalPrefix_(intervalPrefix)
{
auto tControl = dict.getVal<word>("timeControl");
if(tControl == "timeStep")
@ -162,3 +164,32 @@ pFlow::baseTimeControl::iInterval() const
fatalExit;
return 0;
}
bool
pFlow::baseTimeControl::write(dictionary& dict) const
{
if(isTimeStep_)
{
dict.add("timeControl", "timeStep");
}
else
{
dict.add("timeControl", "runTime");
}
word intervalWord = intervalPrefix_.size()==0? word("interval"): intervalPrefix_+"Interval";
if(!isTimeStep_)
{
dict.add(intervalWord,rRange_.stride());
dict.add("startTime",rRange_.begin());
dict.add("endTime", rRange_.end());
}
else
{
dict.add(intervalWord,iRange_.stride());
dict.add("startTime",iRange_.begin());
dict.add("endTime", iRange_.end());
}
return true;
}

View File

@ -1,20 +1,20 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it
under the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#ifndef __baseTimeControl_hpp__
@ -23,50 +23,51 @@ Licence:
#include "dictionary.hpp"
#include "ranges.hpp"
namespace pFlow
{
class baseTimeControl
{
private:
bool isTimeStep_;
bool isTimeStep_;
int32StridedRagne iRange_;
realStridedRange rRange_;
const word intervalPrefix_;
public:
baseTimeControl(
const dictionary& dict,
const word& intervalPrefix = "",
real defStartTime = 0.0);
inline
bool isTimeStep()const
const dictionary& dict,
const word& intervalPrefix = "",
real defStartTime = 0.0
);
inline bool isTimeStep() const
{
return isTimeStep_;
}
bool timeEvent(uint32 iter, real t, real dt)const;
bool timeEvent(uint32 iter, real t, real dt) const;
bool isInRange(uint32 iter, real t, real dt)const;
bool isInRange(uint32 iter, real t, real dt) const;
real startTime()const;
real startTime() const;
real endTime()const;
real endTime() const;
real rInterval()const;
real rInterval() const;
int32 startIter()const;
int32 startIter() const;
int32 endIter()const;
int32 endIter() const;
int32 iInterval()const;
int32 iInterval() const;
bool write(dictionary& dict) const;
};
}

View File

@ -21,11 +21,13 @@ Licence:
#include "cylinder.hpp"
#include "zAxis.hpp"
#include "streams.hpp"
FUNCTION_H
bool pFlow::cylinder::calculateParams()
{
WARNING<<"Use of cylinder requires modifications to zAxis"<<END_WARNING;
auto p1p2 = p2_ - p1_;
if( p1p2.length() > smallValue )

View File

@ -50,17 +50,17 @@ struct triple
T y_;
T z_;
/// Type info for triple
TripleTypeInfoNV(T);
//// Constructors
/// Initilize to zero
INLINE_FUNCTION_HD
triple()
: x_(),
y_(),
z_()
{
}
triple() = default;
INLINE_FUNCTION_HD
~triple() = default;
/// Construct from x, y, z
INLINE_FUNCTION_HD