MPI-boundaries for processor
This commit is contained in:
Hamidreza Norouzi
parent
94fcc3d01b
commit
5f90605a41
|
|
@ -0,0 +1,108 @@
|
|||
/*------------------------------- phasicFlow ---------------------------------
|
||||
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
|
||||
------------------------------------------------------------------------------
|
||||
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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
|
||||
#include "processorBoundaryContactSearch.hpp"
|
||||
#include "contactSearch.hpp"
|
||||
#include "particles.hpp"
|
||||
//#include "pointStructure.hpp"
|
||||
//#include "geometry.hpp"
|
||||
|
||||
|
||||
void pFlow::processorBoundaryContactSearch::setSearchBox()
|
||||
{
|
||||
|
||||
auto l = boundary().neighborLength();
|
||||
auto n = boundary().boundaryPlane().normal();
|
||||
auto pp1 = boundary().boundaryPlane().parallelPlane(l);
|
||||
auto pp2 = boundary().boundaryPlane().parallelPlane(-l);
|
||||
|
||||
realx3 minP1 = min(min(min(pp1.p1(), pp1.p2()), pp1.p3()), pp1.p4());
|
||||
realx3 maxP1 = max(max(max(pp1.p1(), pp1.p2()), pp1.p3()), pp1.p4());
|
||||
|
||||
realx3 minP2 = min(min(min(pp2.p1(), pp2.p2()), pp2.p3()), pp2.p4());
|
||||
realx3 maxP2 = max(max(max(pp2.p1(), pp2.p2()), pp2.p3()), pp2.p4());
|
||||
|
||||
auto minP = min(minP1, minP2) - l*(realx3(1.0)-abs(n));
|
||||
auto maxP = max(maxP1, maxP2) + l*(realx3(1.0)-abs(n));
|
||||
|
||||
searchBox_={minP, maxP};
|
||||
}
|
||||
|
||||
pFlow::processorBoundaryContactSearch::processorBoundaryContactSearch(
|
||||
const dictionary &dict,
|
||||
const boundaryBase &boundary,
|
||||
const contactSearch &cSearch)
|
||||
:
|
||||
boundaryContactSearch(dict, boundary, cSearch),
|
||||
diameter_(cSearch.Particles().boundingSphere()),
|
||||
masterSearch_(this->isBoundaryMaster())
|
||||
{
|
||||
|
||||
if(masterSearch_)
|
||||
{
|
||||
setSearchBox();
|
||||
|
||||
real minD;
|
||||
real maxD;
|
||||
cSearch.Particles().boundingSphereMinMax(minD, maxD);
|
||||
|
||||
ppContactSearch_ = makeUnique<twoPartContactSearch>(
|
||||
searchBox_,
|
||||
maxD);
|
||||
}
|
||||
else
|
||||
{
|
||||
searchBox_={{0,0,0},{0,0,0}};
|
||||
}
|
||||
}
|
||||
|
||||
bool pFlow::processorBoundaryContactSearch::broadSearch
|
||||
(
|
||||
uint32 iter,
|
||||
real t,
|
||||
real dt,
|
||||
csPairContainerType &ppPairs,
|
||||
csPairContainerType &pwPairs,
|
||||
bool force
|
||||
)
|
||||
{
|
||||
if(masterSearch_)
|
||||
{
|
||||
/*const auto thisPoints = boundary().thisPoints();
|
||||
const auto& neighborProcPoints = boundary().neighborProcPoints();
|
||||
const auto& bDiams = diameter_.BoundaryField(thisBoundaryIndex());
|
||||
const auto thisDiams = bDiams.thisField();
|
||||
const auto& neighborProcDiams = bDiams.neighborProcField();
|
||||
|
||||
ppContactSearch_().broadSearchPP(
|
||||
ppPairs,
|
||||
thisPoints,
|
||||
thisDiams,
|
||||
neighborProcPoints,
|
||||
neighborProcDiams);
|
||||
|
||||
pOutput<<"ppPairs size in boundary"<< ppPairs.size()<<endl; */
|
||||
return true;
|
||||
|
||||
}else
|
||||
{
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,74 @@
|
|||
/*------------------------------- phasicFlow ---------------------------------
|
||||
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
|
||||
------------------------------------------------------------------------------
|
||||
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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
#ifndef __processorBoundaryContactSearch_hpp__
|
||||
#define __processorBoundaryContactSearch_hpp__
|
||||
|
||||
#include "boundaryContactSearch.hpp"
|
||||
#include "pointFields.hpp"
|
||||
#include "twoPartContactSearch.hpp"
|
||||
|
||||
namespace pFlow
|
||||
{
|
||||
|
||||
class processorBoundaryContactSearch : public boundaryContactSearch
|
||||
{
|
||||
private:
|
||||
|
||||
box searchBox_;
|
||||
|
||||
uniquePtr<twoPartContactSearch> ppContactSearch_ = nullptr;
|
||||
|
||||
const realPointField_D& diameter_;
|
||||
|
||||
bool masterSearch_;
|
||||
|
||||
void setSearchBox();
|
||||
|
||||
public:
|
||||
|
||||
TypeInfo("boundaryContactSearch<MPI,processor>")
|
||||
|
||||
processorBoundaryContactSearch(
|
||||
const dictionary& dict,
|
||||
const boundaryBase& boundary,
|
||||
const contactSearch& cSearch
|
||||
);
|
||||
|
||||
~processorBoundaryContactSearch() override = default;
|
||||
|
||||
add_vCtor(
|
||||
boundaryContactSearch,
|
||||
processorBoundaryContactSearch,
|
||||
boundaryBase
|
||||
);
|
||||
|
||||
bool broadSearch(
|
||||
uint32 iter,
|
||||
real t,
|
||||
real dt,
|
||||
csPairContainerType& ppPairs,
|
||||
csPairContainerType& pwPairs,
|
||||
bool force = false
|
||||
) override;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif //__processorBoundaryContactSearch_hpp__
|
||||
|
|
@ -0,0 +1,160 @@
|
|||
|
||||
#include "twoPartContactSearch.hpp"
|
||||
#include "twoPartContactSearchKernels.hpp"
|
||||
#include "phasicFlowKokkos.hpp"
|
||||
#include "streams.hpp"
|
||||
|
||||
void pFlow::twoPartContactSearch::checkAllocateNext(uint32 n)
|
||||
{
|
||||
if( nextCapacity_ < n)
|
||||
{
|
||||
nextCapacity_ = n;
|
||||
reallocNoInit(next_, n);
|
||||
}
|
||||
}
|
||||
|
||||
void pFlow::twoPartContactSearch::nullifyHead()
|
||||
{
|
||||
fill(head_, static_cast<uint32>(-1));
|
||||
}
|
||||
|
||||
void pFlow::twoPartContactSearch::nullifyNext(uint32 n)
|
||||
{
|
||||
fill(next_, 0u, n, static_cast<uint32>(-1));
|
||||
}
|
||||
|
||||
void pFlow::twoPartContactSearch::buildList(
|
||||
const deviceScatteredFieldAccess<realx3> &points)
|
||||
{
|
||||
if(points.empty())return;
|
||||
uint32 n = points.size();
|
||||
checkAllocateNext(n);
|
||||
nullifyNext(n);
|
||||
nullifyHead();
|
||||
|
||||
pFlow::twoPartContactSearchKernels::buildNextHead(
|
||||
points,
|
||||
searchCells_,
|
||||
head_,
|
||||
next_
|
||||
);
|
||||
}
|
||||
|
||||
pFlow::twoPartContactSearch::twoPartContactSearch
|
||||
(
|
||||
const box &domain,
|
||||
real cellSize,
|
||||
real sizeRatio
|
||||
)
|
||||
:
|
||||
searchCells_(domain, cellSize),
|
||||
head_("periodic:head",searchCells_.nx(), searchCells_.ny(), searchCells_.nz()),
|
||||
sizeRatio_(sizeRatio)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
bool pFlow::twoPartContactSearch::broadSearchPP
|
||||
(
|
||||
csPairContainerType &ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3> &points1,
|
||||
const deviceScatteredFieldAccess<real>& diams1,
|
||||
const deviceScatteredFieldAccess<realx3> &points2,
|
||||
const deviceScatteredFieldAccess<real>& diams2,
|
||||
const realx3& transferVec
|
||||
)
|
||||
{
|
||||
|
||||
buildList(points1);
|
||||
|
||||
uint32 nNotInserted = 1;
|
||||
|
||||
// loop until the container size fits the numebr of contact pairs
|
||||
while (nNotInserted > 0)
|
||||
{
|
||||
|
||||
nNotInserted = pFlow::twoPartContactSearchKernels::broadSearchPP
|
||||
(
|
||||
ppPairs,
|
||||
points1,
|
||||
diams1,
|
||||
points2,
|
||||
diams2,
|
||||
transferVec,
|
||||
head_,
|
||||
next_,
|
||||
searchCells_,
|
||||
sizeRatio_
|
||||
);
|
||||
|
||||
|
||||
if(nNotInserted)
|
||||
{
|
||||
// - resize the container
|
||||
// note that getFull now shows the number of failed insertions.
|
||||
uint32 len = max(nNotInserted,100u) ;
|
||||
|
||||
auto oldCap = ppPairs.capacity();
|
||||
|
||||
ppPairs.increaseCapacityBy(len);
|
||||
|
||||
INFORMATION<< "Particle-particle contact pair container capacity increased from "<<
|
||||
oldCap << " to "<<ppPairs.capacity()<<" in peiodicBoundaryContactSearch."<<END_INFO;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
bool pFlow::twoPartContactSearch::broadSearchPP
|
||||
(
|
||||
csPairContainerType &ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3> &points1,
|
||||
const deviceScatteredFieldAccess<real> &diams1,
|
||||
const realx3Vector_D& points2,
|
||||
const realVector_D& diams2
|
||||
)
|
||||
{
|
||||
buildList(points1);
|
||||
|
||||
uint32 nNotInserted = 1;
|
||||
|
||||
// loop until the container size fits the numebr of contact pairs
|
||||
while (nNotInserted > 0)
|
||||
{
|
||||
|
||||
nNotInserted = pFlow::twoPartContactSearchKernels::broadSearchPP
|
||||
(
|
||||
ppPairs,
|
||||
points1,
|
||||
diams1,
|
||||
points2,
|
||||
diams2,
|
||||
head_,
|
||||
next_,
|
||||
searchCells_,
|
||||
sizeRatio_
|
||||
);
|
||||
|
||||
|
||||
if(nNotInserted)
|
||||
{
|
||||
// - resize the container
|
||||
// note that getFull now shows the number of failed insertions.
|
||||
uint32 len = max(nNotInserted,100u) ;
|
||||
|
||||
auto oldCap = ppPairs.capacity();
|
||||
|
||||
ppPairs.increaseCapacityBy(len);
|
||||
|
||||
INFORMATION<< "Particle-particle contact pair container capacity increased from "<<
|
||||
oldCap << " to "<<ppPairs.capacity()<<" in peiodicBoundaryContactSearch."<<END_INFO;
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
|
@ -0,0 +1,103 @@
|
|||
/*------------------------------- phasicFlow ---------------------------------
|
||||
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
|
||||
------------------------------------------------------------------------------
|
||||
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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef __twoPartContactSearch_hpp__
|
||||
#define __twoPartContactSearch_hpp__
|
||||
|
||||
#include "contactSearchGlobals.hpp"
|
||||
#include "scatteredFieldAccess.hpp"
|
||||
#include "cells.hpp"
|
||||
#include "VectorSingles.hpp"
|
||||
|
||||
namespace pFlow
|
||||
{
|
||||
|
||||
class twoPartContactSearch
|
||||
{
|
||||
public:
|
||||
using HeadType = deviceViewType3D<uint32>;
|
||||
|
||||
using NextType = deviceViewType1D<uint32>;
|
||||
|
||||
private:
|
||||
|
||||
cells searchCells_;
|
||||
|
||||
HeadType head_{ "periodic::head", 1, 1, 1 };
|
||||
|
||||
NextType next_{ "periodic::next", 1 };
|
||||
|
||||
real sizeRatio_ = 1.0;
|
||||
|
||||
uint32 nextCapacity_ = 0;
|
||||
|
||||
void checkAllocateNext(uint32 n);
|
||||
|
||||
void nullifyHead();
|
||||
|
||||
void nullifyNext(uint32 n);
|
||||
|
||||
void buildList(
|
||||
const deviceScatteredFieldAccess<realx3> &points);
|
||||
|
||||
public:
|
||||
twoPartContactSearch(
|
||||
const box &domain,
|
||||
real cellSize,
|
||||
real sizeRatio = 1.0);
|
||||
|
||||
/// @brief Perform a broad-search for spheres in two adjacent regions.
|
||||
/// Region 1 is considered as the master (primary) region and region 2 as slave
|
||||
/// @param ppPairs pairs container which holds i and j
|
||||
/// @param points1 point positions in region 1
|
||||
/// @param diams1 diameter of spheres in region 1
|
||||
/// @param points2 point positions in region 2
|
||||
/// @param diams2 diameter of spheres in region 2
|
||||
/// @param transferVec a vector to transfer points from region 2 to region 1
|
||||
/// @return true if it is successful
|
||||
bool broadSearchPP(
|
||||
csPairContainerType &ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3> &points1,
|
||||
const deviceScatteredFieldAccess<real> &diams1,
|
||||
const deviceScatteredFieldAccess<realx3> &points2,
|
||||
const deviceScatteredFieldAccess<real> &diams2,
|
||||
const realx3 &transferVec);
|
||||
|
||||
bool broadSearchPP(
|
||||
csPairContainerType &ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3> &points1,
|
||||
const deviceScatteredFieldAccess<real> &diams1,
|
||||
const realx3Vector_D& points2,
|
||||
const realVector_D& diams2);
|
||||
|
||||
const auto& searchCells()const
|
||||
{
|
||||
return searchCells_;
|
||||
}
|
||||
|
||||
real sizeRatio()const
|
||||
{
|
||||
return sizeRatio_;
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif //__twoPartContactSearch_hpp__
|
||||
|
|
@ -0,0 +1,188 @@
|
|||
#include "twoPartContactSearchKernels.hpp"
|
||||
|
||||
INLINE_FUNCTION_HD
|
||||
bool
|
||||
sphereSphereCheckB(
|
||||
const pFlow::realx3& p1,
|
||||
const pFlow::realx3 p2,
|
||||
pFlow::real d1,
|
||||
pFlow::real d2
|
||||
)
|
||||
{
|
||||
return pFlow::length(p2 - p1) < 0.5 * (d2 + d1);
|
||||
}
|
||||
|
||||
void
|
||||
pFlow::twoPartContactSearchKernels::buildNextHead(
|
||||
const deviceScatteredFieldAccess<realx3>& points,
|
||||
const cells& searchCells,
|
||||
deviceViewType3D<uint32>& head,
|
||||
deviceViewType1D<uint32>& next
|
||||
)
|
||||
{
|
||||
if (points.empty())
|
||||
return;
|
||||
|
||||
uint32 n = points.size();
|
||||
|
||||
Kokkos::parallel_for(
|
||||
"pFlow::ppwBndryContactSearch::buildList",
|
||||
deviceRPolicyStatic(0, n),
|
||||
LAMBDA_HD(uint32 i) {
|
||||
int32x3 ind;
|
||||
if (searchCells.pointIndexInDomain(points[i], ind))
|
||||
{
|
||||
// discards points out of searchCell
|
||||
uint32 old =
|
||||
Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
|
||||
next[i] = old;
|
||||
}
|
||||
}
|
||||
);
|
||||
Kokkos::fence();
|
||||
}
|
||||
|
||||
pFlow::uint32
|
||||
pFlow::twoPartContactSearchKernels::broadSearchPP(
|
||||
csPairContainerType& ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3>& points,
|
||||
const deviceScatteredFieldAccess<real>& diams,
|
||||
const deviceScatteredFieldAccess<realx3>& mirrorPoints,
|
||||
const deviceScatteredFieldAccess<real>& mirrorDiams,
|
||||
const realx3& transferVec,
|
||||
const deviceViewType3D<uint32>& head,
|
||||
const deviceViewType1D<uint32>& next,
|
||||
const cells& searchCells,
|
||||
const real sizeRatio
|
||||
)
|
||||
{
|
||||
if (points.empty())
|
||||
return 0;
|
||||
if (mirrorPoints.empty())
|
||||
return 0;
|
||||
|
||||
auto nMirror = mirrorPoints.size();
|
||||
|
||||
uint32 getFull = 0;
|
||||
|
||||
Kokkos::parallel_reduce(
|
||||
"pFlow::twoPartContactSearchKernels::broadSearchPP",
|
||||
deviceRPolicyStatic(0, nMirror),
|
||||
LAMBDA_HD(const uint32 mrrI, uint32& getFullUpdate) {
|
||||
realx3 p_m = mirrorPoints(mrrI) + transferVec;
|
||||
|
||||
int32x3 ind_m;
|
||||
if (!searchCells.pointIndexInDomain(p_m, ind_m))
|
||||
return;
|
||||
|
||||
real d_m = sizeRatio * mirrorDiams[mrrI];
|
||||
|
||||
for (int ii = -1; ii < 2; ii++)
|
||||
{
|
||||
for (int jj = -1; jj < 2; jj++)
|
||||
{
|
||||
for (int kk = -1; kk < 2; kk++)
|
||||
{
|
||||
auto ind = ind_m + int32x3{ ii, jj, kk };
|
||||
|
||||
if (!searchCells.inCellRange(ind))
|
||||
continue;
|
||||
|
||||
uint32 thisI = head(ind.x(), ind.y(), ind.z());
|
||||
while (thisI != -1)
|
||||
{
|
||||
auto d_n = sizeRatio * diams[thisI];
|
||||
|
||||
// first item is for this boundary and second itme,
|
||||
// for mirror
|
||||
if(sphereSphereCheckB(p_m, points[thisI], d_m, d_n)&&
|
||||
ppPairs.insert(thisI,mrrI) == -1)
|
||||
{
|
||||
getFullUpdate++;
|
||||
}
|
||||
|
||||
thisI = next(thisI);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
},
|
||||
getFull
|
||||
);
|
||||
|
||||
return getFull;
|
||||
}
|
||||
|
||||
pFlow::uint32
|
||||
pFlow::twoPartContactSearchKernels::broadSearchPP(
|
||||
csPairContainerType& ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3>& points1,
|
||||
const deviceScatteredFieldAccess<real>& diams1,
|
||||
const realx3Vector_D& points2,
|
||||
const realVector_D& diams2,
|
||||
const deviceViewType3D<uint32>& head,
|
||||
const deviceViewType1D<uint32>& next,
|
||||
const cells& searchCells,
|
||||
real sizeRatio
|
||||
)
|
||||
{
|
||||
if (points1.empty())
|
||||
return 0;
|
||||
if (points2.empty())
|
||||
return 0;
|
||||
|
||||
auto nP2 = points2.size();
|
||||
auto points2View = points2.deviceView();
|
||||
auto diams2View = diams2.deviceView();
|
||||
|
||||
uint32 getFull = 0;
|
||||
|
||||
Kokkos::parallel_reduce(
|
||||
"pFlow::twoPartContactSearchKernels::broadSearchPP",
|
||||
deviceRPolicyStatic(0, nP2),
|
||||
LAMBDA_HD(const uint32 i2, uint32& getFullUpdate) {
|
||||
realx3 p_m = points2View(i2);
|
||||
|
||||
int32x3 ind_m;
|
||||
if (!searchCells.pointIndexInDomain(p_m, ind_m))
|
||||
return;
|
||||
|
||||
real d_m = sizeRatio * diams2View[i2];
|
||||
|
||||
for (int ii = -1; ii < 2; ii++)
|
||||
{
|
||||
for (int jj = -1; jj < 2; jj++)
|
||||
{
|
||||
for (int kk = -1; kk < 2; kk++)
|
||||
{
|
||||
auto ind = ind_m + int32x3{ ii, jj, kk };
|
||||
|
||||
if (!searchCells.inCellRange(ind))
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
uint32 i1 = head(ind.x(), ind.y(), ind.z());
|
||||
while (i1 != -1)
|
||||
{
|
||||
auto d_n = sizeRatio * diams1[i1];
|
||||
|
||||
// first item is for this boundary and second itme,
|
||||
// for mirror
|
||||
if(sphereSphereCheckB(p_m, points1[i1], d_m, d_n)&&
|
||||
ppPairs.insert(i1,i2) == -1)
|
||||
{
|
||||
getFullUpdate++;
|
||||
}
|
||||
|
||||
i1 = next(i1);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
},
|
||||
getFull
|
||||
);
|
||||
|
||||
return getFull;
|
||||
}
|
||||
|
|
@ -0,0 +1,49 @@
|
|||
#ifndef __twoPartContactSearchKernels_hpp__
|
||||
#define __twoPartContactSearchKernels_hpp__
|
||||
|
||||
#include "contactSearchGlobals.hpp"
|
||||
#include "cells.hpp"
|
||||
#include "contactSearchFunctions.hpp"
|
||||
#include "scatteredFieldAccess.hpp"
|
||||
#include "VectorSingles.hpp"
|
||||
|
||||
namespace pFlow::twoPartContactSearchKernels
|
||||
{
|
||||
|
||||
void buildNextHead(
|
||||
const deviceScatteredFieldAccess<realx3> &points,
|
||||
const cells &searchCells,
|
||||
deviceViewType3D<uint32> &head,
|
||||
deviceViewType1D<uint32> &next );
|
||||
|
||||
|
||||
uint32 broadSearchPP
|
||||
(
|
||||
csPairContainerType &ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3> &points,
|
||||
const deviceScatteredFieldAccess<real> &diams,
|
||||
const deviceScatteredFieldAccess<realx3> &mirrorPoints,
|
||||
const deviceScatteredFieldAccess<real> &mirrorDiams,
|
||||
const realx3 &transferVec,
|
||||
const deviceViewType3D<uint32> &head,
|
||||
const deviceViewType1D<uint32> &next,
|
||||
const cells &searchCells,
|
||||
real sizeRatio
|
||||
);
|
||||
|
||||
uint32
|
||||
broadSearchPP(
|
||||
csPairContainerType& ppPairs,
|
||||
const deviceScatteredFieldAccess<realx3>& points1,
|
||||
const deviceScatteredFieldAccess<real>& diams1,
|
||||
const realx3Vector_D& points2,
|
||||
const realVector_D& diams2,
|
||||
const deviceViewType3D<uint32>& head,
|
||||
const deviceViewType1D<uint32>& next,
|
||||
const cells& searchCells,
|
||||
real sizeRatio
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
#endif //__twoPartContactSearchKernels_hpp__
|
||||
|
|
@ -0,0 +1,131 @@
|
|||
|
||||
#ifndef __processorBoundarySIKernels_hpp__
|
||||
#define __processorBoundarySIKernels_hpp__
|
||||
|
||||
namespace pFlow::MPI::processorBoundarySIKernels
|
||||
{
|
||||
|
||||
template<typename ContactListType, typename ContactForceModel>
|
||||
inline
|
||||
void sphereSphereInteraction
|
||||
(
|
||||
real dt,
|
||||
const ContactListType& cntctList,
|
||||
const ContactForceModel& forceModel,
|
||||
const deviceScatteredFieldAccess<realx3>& thisPoints,
|
||||
const deviceViewType1D<real>& thisDiam,
|
||||
const deviceViewType1D<uint32>& thisPropId,
|
||||
const deviceViewType1D<realx3>& thisVel,
|
||||
const deviceViewType1D<realx3>& thisRVel,
|
||||
const deviceViewType1D<realx3>& thisCForce,
|
||||
const deviceViewType1D<realx3>& thisCTorque,
|
||||
const deviceViewType1D<realx3>& neighborPoints,
|
||||
const deviceViewType1D<real>& neighborDiam,
|
||||
const deviceViewType1D<uint32>& neighborPropId,
|
||||
const deviceViewType1D<realx3>& neighborVel,
|
||||
const deviceViewType1D<realx3>& neighborRVel,
|
||||
const deviceViewType1D<realx3>& neighborCForce,
|
||||
const deviceViewType1D<realx3>& neighborCTorque
|
||||
)
|
||||
{
|
||||
|
||||
using ValueType = typename ContactListType::ValueType;
|
||||
uint32 ss = cntctList.size();
|
||||
if(ss == 0u)return;
|
||||
|
||||
uint32 lastItem = cntctList.loopCount();
|
||||
|
||||
Kokkos::parallel_for(
|
||||
"pFlow::MPI::processorBoundarySIKernels::sphereSphereInteraction",
|
||||
deviceRPolicyDynamic(0,lastItem),
|
||||
LAMBDA_HD(uint32 n)
|
||||
{
|
||||
|
||||
if(!cntctList.isValid(n))return;
|
||||
|
||||
auto [i,j] = cntctList.getPair(n);
|
||||
uint32 ind_i = thisPoints.index(i);
|
||||
uint32 ind_j = j;
|
||||
|
||||
real Ri = 0.5*thisDiam[ind_i];
|
||||
real Rj = 0.5*neighborDiam[ind_j];
|
||||
realx3 xi = thisPoints.field()[ind_i];
|
||||
realx3 xj = neighborPoints[ind_j];
|
||||
|
||||
real dist = length(xj-xi);
|
||||
real ovrlp = (Ri+Rj) - dist;
|
||||
|
||||
if( ovrlp >0.0 )
|
||||
{
|
||||
auto Nij = (xj-xi)/max(dist,smallValue);
|
||||
auto wi = thisRVel[ind_i];
|
||||
auto wj = neighborRVel[ind_j];
|
||||
auto Vr = thisVel[ind_i] - neighborVel[ind_j] + cross((Ri*wi+Rj*wj), Nij);
|
||||
|
||||
auto history = cntctList.getValue(n);
|
||||
|
||||
int32 propId_i = thisPropId[ind_i];
|
||||
int32 propId_j = neighborPropId[ind_j];
|
||||
|
||||
realx3 FCn, FCt, Mri, Mrj, Mij, Mji;
|
||||
|
||||
// calculates contact force
|
||||
forceModel.contactForce(
|
||||
dt, i, j,
|
||||
propId_i, propId_j,
|
||||
Ri, Rj,
|
||||
ovrlp,
|
||||
Vr, Nij,
|
||||
history,
|
||||
FCn, FCt);
|
||||
|
||||
forceModel.rollingFriction(
|
||||
dt, i, j,
|
||||
propId_i, propId_j,
|
||||
Ri, Rj,
|
||||
wi, wj,
|
||||
Nij,
|
||||
FCn,
|
||||
Mri, Mrj);
|
||||
|
||||
auto M = cross(Nij,FCt);
|
||||
Mij = Ri*M+Mri;
|
||||
Mji = Rj*M+Mrj;
|
||||
|
||||
auto FC = FCn + FCt;
|
||||
|
||||
|
||||
Kokkos::atomic_add(&thisCForce[ind_i].x_,FC.x_);
|
||||
Kokkos::atomic_add(&thisCForce[ind_i].y_,FC.y_);
|
||||
Kokkos::atomic_add(&thisCForce[ind_i].z_,FC.z_);
|
||||
|
||||
Kokkos::atomic_add(&neighborCForce[ind_j].x_,-FC.x_);
|
||||
Kokkos::atomic_add(&neighborCForce[ind_j].y_,-FC.y_);
|
||||
Kokkos::atomic_add(&neighborCForce[ind_j].z_,-FC.z_);
|
||||
|
||||
Kokkos::atomic_add(&thisCTorque[ind_i].x_, Mij.x_);
|
||||
Kokkos::atomic_add(&thisCTorque[ind_i].y_, Mij.y_);
|
||||
Kokkos::atomic_add(&thisCTorque[ind_i].z_, Mij.z_);
|
||||
|
||||
Kokkos::atomic_add(&neighborCTorque[ind_j].x_, Mji.x_);
|
||||
Kokkos::atomic_add(&neighborCTorque[ind_j].y_, Mji.y_);
|
||||
Kokkos::atomic_add(&neighborCTorque[ind_j].z_, Mji.z_);
|
||||
|
||||
|
||||
cntctList.setValue(n,history);
|
||||
|
||||
}
|
||||
else
|
||||
{
|
||||
cntctList.setValue(n, ValueType());
|
||||
}
|
||||
|
||||
});
|
||||
Kokkos::fence();
|
||||
}
|
||||
|
||||
|
||||
} //pFlow::MPI::processorBoundarySIKernels
|
||||
|
||||
|
||||
#endif //__processorBoundarySIKernels_hpp__
|
||||
|
|
@ -0,0 +1,73 @@
|
|||
/*------------------------------- phasicFlow ---------------------------------
|
||||
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
|
||||
------------------------------------------------------------------------------
|
||||
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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
|
||||
#include "processorBoundarySIKernels.hpp"
|
||||
|
||||
template <typename cFM, typename gMM>
|
||||
pFlow::MPI::processorBoundarySphereInteraction<cFM, gMM>::processorBoundarySphereInteraction(
|
||||
const boundaryBase &boundary,
|
||||
const sphereParticles &sphPrtcls,
|
||||
const GeometryMotionModel &geomMotion)
|
||||
:
|
||||
boundarySphereInteraction<cFM,gMM>(
|
||||
boundary,
|
||||
sphPrtcls,
|
||||
geomMotion
|
||||
),
|
||||
masterInteraction_(boundary.isBoundaryMaster())
|
||||
{
|
||||
pOutput<<"Processor boundayrCondition for "<< boundary.name()<<endl;
|
||||
}
|
||||
|
||||
template <typename cFM, typename gMM>
|
||||
bool pFlow::MPI::processorBoundarySphereInteraction<cFM, gMM>::sphereSphereInteraction
|
||||
(
|
||||
real dt,
|
||||
const ContactForceModel &cfModel
|
||||
)
|
||||
{
|
||||
if(!masterInteraction_) return true;
|
||||
|
||||
const auto & sphPar = this->sphParticles();
|
||||
uint32 thisIndex = this->boundary().thisBoundaryIndex();
|
||||
const auto& a = sphPar.diameter().BoundaryField(thisIndex).neighborProcField().deviceViewAll();
|
||||
|
||||
/*pFlow::MPI::processorBoundarySIKernels::sphereSphereInteraction(
|
||||
dt,
|
||||
this->ppPairs(),
|
||||
cfModel,
|
||||
this->boundary().thisPoints(),
|
||||
sphPar.diameter().deviceViewAll(),
|
||||
sphPar.propertyId().deviceViewAll(),
|
||||
sphPar.velocity().deviceViewAll(),
|
||||
sphPar.rVelocity().deviceViewAll(),
|
||||
sphPar.contactForce().deviceViewAll(),
|
||||
sphPar.contactTorque().deviceViewAll(),
|
||||
this->boundary().neighborProcPoints().deviceViewAll(),
|
||||
sphPar.diameter().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
|
||||
sphPar.propertyId().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
|
||||
sphPar.velocity().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
|
||||
sphPar.rVelocity().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
|
||||
sphPar.contactForce().BoundaryField(thisIndex).neighborProcField().deviceViewAll(),
|
||||
sphPar.contactTorque().BoundaryField(thisIndex).neighborProcField().deviceViewAll()
|
||||
);*/
|
||||
|
||||
return true;
|
||||
}
|
||||
|
|
@ -0,0 +1,90 @@
|
|||
/*------------------------------- phasicFlow ---------------------------------
|
||||
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
|
||||
------------------------------------------------------------------------------
|
||||
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.
|
||||
|
||||
-----------------------------------------------------------------------------*/
|
||||
#ifndef __processorBoundarySphereInteraction_hpp__
|
||||
#define __processorBoundarySphereInteraction_hpp__
|
||||
|
||||
#include "boundarySphereInteraction.hpp"
|
||||
|
||||
namespace pFlow::MPI
|
||||
{
|
||||
|
||||
template<typename contactForceModel,typename geometryMotionModel>
|
||||
class processorBoundarySphereInteraction
|
||||
:
|
||||
public boundarySphereInteraction<contactForceModel, geometryMotionModel>
|
||||
{
|
||||
public:
|
||||
|
||||
using PBSInteractionType =
|
||||
processorBoundarySphereInteraction<contactForceModel,geometryMotionModel>;
|
||||
|
||||
using BSInteractionType =
|
||||
boundarySphereInteraction<contactForceModel, geometryMotionModel>;
|
||||
|
||||
using GeometryMotionModel = typename BSInteractionType::GeometryMotionModel;
|
||||
|
||||
using ContactForceModel = typename BSInteractionType::ContactForceModel;
|
||||
|
||||
using MotionModel = typename geometryMotionModel::MotionModel;
|
||||
|
||||
using ModelStorage = typename ContactForceModel::contactForceStorage;
|
||||
|
||||
using IdType = typename BSInteractionType::IdType;
|
||||
|
||||
using IndexType = typename BSInteractionType::IndexType;
|
||||
|
||||
using ContactListType = typename BSInteractionType::ContactListType;
|
||||
|
||||
private:
|
||||
|
||||
bool masterInteraction_;
|
||||
|
||||
public:
|
||||
|
||||
TypeInfoTemplate22("boundarySphereInteraction", "processor",ContactForceModel, MotionModel);
|
||||
|
||||
|
||||
processorBoundarySphereInteraction(
|
||||
const boundaryBase& boundary,
|
||||
const sphereParticles& sphPrtcls,
|
||||
const GeometryMotionModel& geomMotion
|
||||
);
|
||||
|
||||
add_vCtor
|
||||
(
|
||||
BSInteractionType,
|
||||
PBSInteractionType,
|
||||
boundaryBase
|
||||
);
|
||||
|
||||
~processorBoundarySphereInteraction()override = default;
|
||||
|
||||
bool sphereSphereInteraction(
|
||||
real dt,
|
||||
const ContactForceModel& cfModel)override;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#include "processorBoundarySphereInteraction.cpp"
|
||||
|
||||
|
||||
#endif //__processorBoundarySphereInteraction_hpp__
|
||||
|
|
@ -0,0 +1,17 @@
|
|||
|
||||
#include "processorBoundarySphereInteraction.hpp"
|
||||
#include "geometryMotions.hpp"
|
||||
#include "contactForceModels.hpp"
|
||||
|
||||
|
||||
template class pFlow::MPI::processorBoundarySphereInteraction
|
||||
<
|
||||
pFlow::cfModels::limitedNonLinearNormalRolling,
|
||||
pFlow::rotationAxisMotionGeometry
|
||||
>;
|
||||
|
||||
template class pFlow::MPI::processorBoundarySphereInteraction
|
||||
<
|
||||
pFlow::cfModels::nonLimitedNonLinearNormalRolling,
|
||||
pFlow::rotationAxisMotionGeometry
|
||||
>;
|
||||
Loading…
Reference in New Issue