From e10ee2b6b5f10769085bac793bba1dd31f4c5331 Mon Sep 17 00:00:00 2001 From: Hamidreza Norouzi Date: Thu, 1 Feb 2024 12:32:15 -0800 Subject: [PATCH] periodic boundary condition is added to pointStructure --- src/phasicFlow/CMakeLists.txt | 3 + src/phasicFlow/Kokkos/KokkosTypes.hpp | 21 ++- .../containers/VectorHD/VectorSingle.cpp | 9 +- .../containers/VectorHD/VectorSingle.hpp | 6 +- .../pointField/periodicBoundaryField.cpp | 29 ++++ .../pointField/periodicBoundaryField.hpp | 83 ++++++++++ .../containers/pointField/pointFields.cpp | 34 ++-- .../boundaries/boundaryBase/boundaryBase.cpp | 145 ++++++++++++++++-- .../boundaries/boundaryBase/boundaryBase.hpp | 28 +++- .../boundaryBase/boundaryBaseKernels.cpp | 90 +++++++++++ .../boundaryBase/boundaryBaseKernels.hpp | 46 ++++++ .../boundaries/boundaryExit/boundaryExit.cpp | 86 +---------- .../boundaries/boundaryList.cpp | 12 +- .../boundaryPeriodic/boundaryPeriodic.cpp | 125 +++++++++++++++ .../boundaryPeriodic/boundaryPeriodic.hpp | 76 +++++++++ .../domain/regularSimulationDomain.cpp | 24 ++- .../domain/simulationDomain.cpp | 15 +- .../domain/simulationDomain.hpp | 5 + .../pointStructure/internalPoints.cpp | 49 ++++++ .../pointStructure/internalPoints.hpp | 27 ++++ .../pointStructure/internalPointsKernels.cpp | 44 ++++++ .../pointStructure/internalPointsKernels.hpp | 31 ++++ .../pointStructure/pointFlag.hpp | 41 ++++- .../pointStructure/pointFlagKernels.hpp | 78 +++++++++- .../pointStructure/pointStructure.hpp | 16 +- 25 files changed, 988 insertions(+), 135 deletions(-) create mode 100644 src/phasicFlow/containers/pointField/periodicBoundaryField.cpp create mode 100644 src/phasicFlow/containers/pointField/periodicBoundaryField.hpp create mode 100644 src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp create mode 100644 src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.hpp create mode 100644 src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp create mode 100644 src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp create mode 100644 src/phasicFlow/structuredData/pointStructure/internalPointsKernels.cpp create mode 100644 src/phasicFlow/structuredData/pointStructure/internalPointsKernels.hpp diff --git a/src/phasicFlow/CMakeLists.txt b/src/phasicFlow/CMakeLists.txt index 2afe777f..29ec8261 100644 --- a/src/phasicFlow/CMakeLists.txt +++ b/src/phasicFlow/CMakeLists.txt @@ -71,11 +71,14 @@ structuredData/plane/plane.cpp structuredData/domain/domain.cpp structuredData/domain/simulationDomain.cpp structuredData/domain/regularSimulationDomain.cpp +structuredData/pointStructure/internalPointsKernels.cpp structuredData/pointStructure/internalPoints.cpp structuredData/pointStructure/pointStructure.cpp structuredData/boundaries/boundaryBase/boundaryBase.cpp +structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp structuredData/boundaries/boundaryExit/boundaryExit.cpp structuredData/boundaries/boundaryNone/boundaryNone.cpp +structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp structuredData/boundaries/boundaryList.cpp structuredData/pointStructure/selectors/pStructSelector/pStructSelector.cpp structuredData/pointStructure/selectors/selectBox/selectBox.cpp diff --git a/src/phasicFlow/Kokkos/KokkosTypes.hpp b/src/phasicFlow/Kokkos/KokkosTypes.hpp index 3d8cac38..d58df4cb 100644 --- a/src/phasicFlow/Kokkos/KokkosTypes.hpp +++ b/src/phasicFlow/Kokkos/KokkosTypes.hpp @@ -39,13 +39,6 @@ Licence: namespace pFlow { -///********DEP -class DeviceSide{}; -class HostSide{}; - -template -struct selectSide{}; -/*********/// /// Host memory space using HostSpace = Kokkos::HostSpace; @@ -71,6 +64,20 @@ using DefaultHostExecutionSpace = Kokkos::DefaultHostExecutionSpace; using DefaultExecutionSpace = Kokkos::DefaultExecutionSpace; +using deviceRPolicyStatic = + Kokkos::RangePolicy< + Kokkos::DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType >; + + +using hostRPolicyStatic = + Kokkos::RangePolicy< + Kokkos::DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType >; + + /// Pair of two variables template using Pair = Kokkos::pair; diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp index 0fe5d251..451c091a 100644 --- a/src/phasicFlow/containers/VectorHD/VectorSingle.cpp +++ b/src/phasicFlow/containers/VectorHD/VectorSingle.cpp @@ -232,7 +232,14 @@ pFlow::VectorSingle::VectorField()const template INLINE_FUNCTION_H -auto pFlow::VectorSingle::deviceViewAll() const +const auto& pFlow::VectorSingle::deviceViewAll() const +{ + return view_; +} + +template +INLINE_FUNCTION_H +auto& pFlow::VectorSingle::deviceViewAll() { return view_; } diff --git a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp index f0185cb9..7f32e1e3 100644 --- a/src/phasicFlow/containers/VectorHD/VectorSingle.hpp +++ b/src/phasicFlow/containers/VectorHD/VectorSingle.hpp @@ -191,7 +191,11 @@ public: /// Device view range [0,capcity) INLINE_FUNCTION_H - auto deviceViewAll() const; + auto& deviceViewAll(); + + /// Device view range [0,capcity) + INLINE_FUNCTION_H + const auto& deviceViewAll() const; /// Device view range [0, size) INLINE_FUNCTION_H diff --git a/src/phasicFlow/containers/pointField/periodicBoundaryField.cpp b/src/phasicFlow/containers/pointField/periodicBoundaryField.cpp new file mode 100644 index 00000000..16f7f965 --- /dev/null +++ b/src/phasicFlow/containers/pointField/periodicBoundaryField.cpp @@ -0,0 +1,29 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +template + pFlow::periodicBoundaryField::periodicBoundaryField +( + const boundaryBase& boundary, + InternalFieldType& internal +) +: + BoundaryFieldType(boundary, internal) +{} \ No newline at end of file diff --git a/src/phasicFlow/containers/pointField/periodicBoundaryField.hpp b/src/phasicFlow/containers/pointField/periodicBoundaryField.hpp new file mode 100644 index 00000000..9890d55d --- /dev/null +++ b/src/phasicFlow/containers/pointField/periodicBoundaryField.hpp @@ -0,0 +1,83 @@ +/*------------------------------- 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 __periodicBoundaryField_hpp__ +#define __periodicBoundaryField_hpp__ + +#include "boundaryField.hpp" + +namespace pFlow +{ + +template< class T, class MemorySpace = void> +class periodicBoundaryField +: + public boundaryField +{ +public: + + using periodicBoundaryFieldType = periodicBoundaryField; + + using BoundaryFieldType = boundaryField; + + using InternalFieldType = typename BoundaryFieldType::InternalFieldType; + + using memory_space = typename BoundaryFieldType::memory_space; + + using execution_space = typename BoundaryFieldType::execution_space; + + + +public: + + TypeInfo("boundaryField"); + + periodicBoundaryField( + const boundaryBase& boundary, + InternalFieldType& internal); + + + add_vCtor + ( + BoundaryFieldType, + periodicBoundaryFieldType, + boundaryBase + ); + + + bool hearChanges + ( + real t, + real dt, + uint32 iter, + const message& msg, + const anyList& varList + ) override + { + notImplementedFunction; + return false; + } + +}; + +} + +#include "periodicBoundaryField.cpp" + +#endif diff --git a/src/phasicFlow/containers/pointField/pointFields.cpp b/src/phasicFlow/containers/pointField/pointFields.cpp index 8e90019f..09a979b1 100644 --- a/src/phasicFlow/containers/pointField/pointFields.cpp +++ b/src/phasicFlow/containers/pointField/pointFields.cpp @@ -21,69 +21,73 @@ Licence: #include "pointFields.hpp" #include "createBoundaryFields.hpp" +#include "periodicBoundaryField.hpp" +#define createAllBoundary(DataType, MemorySpaceType) \ + template class pFlow::exitBoundaryField; \ + template class pFlow::periodicBoundaryField; template class pFlow::pointField; createBaseBoundary(pFlow::int8, pFlow::HostSpace); -createBoundary(pFlow::int8, pFlow::HostSpace, exit); +createAllBoundary(pFlow::int8, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::int8, void); -createBoundary(pFlow::int8, void, exit); +createAllBoundary(pFlow::int8, void); template class pFlow::pointField; createBaseBoundary(pFlow::uint8, pFlow::HostSpace); -createBoundary(pFlow::uint8, pFlow::HostSpace, exit); +createAllBoundary(pFlow::uint8, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::uint8, void); -createBoundary(pFlow::uint8, void, exit); +createAllBoundary(pFlow::uint8, void); template class pFlow::pointField; createBaseBoundary(pFlow::int32, pFlow::HostSpace); -createBoundary(pFlow::int32, pFlow::HostSpace, exit); +createAllBoundary(pFlow::int32, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::int32, void); -createBoundary(pFlow::int32, void, exit); +createAllBoundary(pFlow::int32, void); template class pFlow::pointField; createBaseBoundary(pFlow::uint32, pFlow::HostSpace); -createBoundary(pFlow::uint32, pFlow::HostSpace, exit); +createAllBoundary(pFlow::uint32, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::uint32, void); -createBoundary(pFlow::uint32, void, exit); +createAllBoundary(pFlow::uint32, void); template class pFlow::pointField; createBaseBoundary(pFlow::real, pFlow::HostSpace); -createBoundary(pFlow::real, pFlow::HostSpace, exit); +createAllBoundary(pFlow::real, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::real, void); -createBoundary(pFlow::real, void, exit); +createAllBoundary(pFlow::real, void); template class pFlow::pointField; createBaseBoundary(pFlow::realx3, pFlow::HostSpace); -createBoundary(pFlow::realx3, pFlow::HostSpace, exit); +createAllBoundary(pFlow::realx3, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::realx3, void); -createBoundary(pFlow::realx3, void, exit); +createAllBoundary(pFlow::realx3, void); template class pFlow::pointField; createBaseBoundary(pFlow::realx4, pFlow::HostSpace); -createBoundary(pFlow::realx4, pFlow::HostSpace, exit); +createAllBoundary(pFlow::realx4, pFlow::HostSpace); template class pFlow::pointField; createBaseBoundary(pFlow::realx4, void); -createBoundary(pFlow::realx4, void, exit); +createAllBoundary(pFlow::realx4, void); template class pFlow::pointField; createBaseBoundary(pFlow::word, pFlow::HostSpace); -createBoundary(pFlow::word, pFlow::HostSpace, exit); \ No newline at end of file +createAllBoundary(pFlow::word, pFlow::HostSpace); \ No newline at end of file diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp index cc11e502..d6821ca9 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.cpp @@ -22,6 +22,8 @@ Licence: #include "dictionary.hpp" #include "internalPoints.hpp" +#include "boundaryBaseKernels.hpp" + /*pFlow::boundaryBase::boundaryBase ( const plane& bplane, @@ -48,21 +50,142 @@ void pFlow::boundaryBase::setNewIndices { auto newSize = newIndices.size(); setSize(static_cast(newSize)); - copy(indexList_.deviceView(), newIndices); + if(newSize>0u) + { + copy(indexList_.deviceView(), newIndices); + } } -pFlow::boundaryBase::boundaryBase( +void pFlow::boundaryBase::appendNewIndices +( + deviceViewType1D newIndices +) +{ + auto s = static_cast(newIndices.size()); + if(s == 0) return; + + uint32 oldS = size(); + uint32 newSize = oldS + s; + + setSize(newSize); + auto appendView = Kokkos::subview( + indexList_.deviceViewAll(), + Kokkos::make_pair(oldS, newSize)); + copy(appendView, newIndices); + + // TODO: notify observers about this change + + // the index list should be sorted + //sort(indexList_.deviceViewAll(), 0, newSize); +} + + + +bool pFlow::boundaryBase::removeIndices +( + uint32 numRemove, + deviceViewType1D removeMask +) +{ + if(removeMask.size() != size()+1 ) + { + fatalErrorInFunction<<"size mismatch"< removeIndices("removeIndices", 1); + deviceViewType1D keepIndices("keepIndices",1); + + pFlow::boundaryBaseKernels::createRemoveKeepIndices + ( + indexList_.deviceView(), + numRemove, + removeMask, + removeIndices, + keepIndices + ); + + if(!internal_.deletePoints(removeIndices)) + { + fatalErrorInFunction<< + "error in deleting points from boundary "<< name()< transferMask, + uint32 transferBoundaryIndex, + realx3 transferVector +) +{ + if(transferMask.size() != size()+1 ) + { + fatalErrorInFunction<<"size mismatch"< transferIndices("transferIndices",1); + deviceViewType1D keepIndices("keepIndices",1); + + pFlow::boundaryBaseKernels::createRemoveKeepIndices + ( + indexList_.deviceView(), + numTransfer, + transferMask, + transferIndices, + keepIndices + ); + + // third, remove the indices from this list + setNewIndices(keepIndices); + + // first, change the flags in the internalPoints + if( !internal_.changePointsFlag( + transferIndices, + transferBoundaryIndex) ) + { + return false; + } + + // second, change the position of points + if(!internal_.changePointsPoisition( + transferIndices, + transferVector)) + { + return false; + } + + // fourth, add the indices to the mirror boundary + internal_.boundary(transferBoundaryIndex). + appendNewIndices(transferIndices); + + return true; +} + +pFlow::boundaryBase::boundaryBase +( const dictionary &dict, const plane &bplane, - internalPoints &internal) - : subscriber(dict.name()), - boundaryPlane_(bplane), - indexList_(groupNames(dict.name(), "indexList")), - neighborLength_(dict.getVal("neighborLength")), - internal_(internal), - mirrorProcessoNo_(dict.getVal("mirrorProcessorNo")), - name_(dict.name()), - type_(dict.getVal("type")) + internalPoints &internal +) +: + subscriber(dict.name()), + boundaryPlane_(bplane), + indexList_(groupNames(dict.name(), "indexList")), + neighborLength_(dict.getVal("neighborLength")), + internal_(internal), + mirrorProcessoNo_(dict.getVal("mirrorProcessorNo")), + name_(dict.name()), + type_(dict.getVal("type")) { } diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp index 4974174d..038114e7 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp @@ -69,6 +69,19 @@ protected: void setNewIndices(deviceViewType1D newIndices); + void appendNewIndices(deviceViewType1D newIndices); + + bool removeIndices( + uint32 numRemove, + deviceViewType1D removeMask); + + bool transferPoints( + uint32 numTransfer, + deviceViewType1D transferMask, + uint32 transferBoundaryIndex, + realx3 transferVector); + + public: TypeInfo("boundaryBase"); @@ -102,12 +115,17 @@ public: (dict, bplane, internal) ); - inline - auto neighborLength()const + virtual + real neighborLength()const { return neighborLength_; } + virtual + realx3 boundaryExtensionLength()const + { + return {0,0,0}; + } const word& type()const { @@ -129,6 +147,7 @@ public: return indexList_.size(); } + virtual const plane& boundaryPlane()const { return boundaryPlane_; @@ -168,11 +187,6 @@ public: return indexList_; } - /*auto& indexList() - { - return indexList_; - }*/ - pointFieldAccessType thisPoints(); virtual diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp new file mode 100644 index 00000000..ff31ba02 --- /dev/null +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.cpp @@ -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. + +-----------------------------------------------------------------------------*/ + +#include "boundaryBaseKernels.hpp" + +void pFlow::boundaryBaseKernels::createRemoveKeepLists +( + uint32 numTotal, + uint32 numRemove, + deviceViewType1D removeMask, + deviceViewType1D& removeList, + deviceViewType1D& keepList +) +{ + uint32 numKeep = numTotal - numRemove; + deviceViewType1D rList("rList",numRemove); + deviceViewType1D kList("kList",numKeep); + + exclusiveScan(removeMask, 0u, numTotal+1, removeMask, 0u); + + Kokkos::parallel_for + ( + "pFlow::boundaryBaseKernels::createRemoveKeepLists", + deviceRPolicyStatic(0, numTotal), + LAMBDA_HD(uint32 i) + { + if(removeMask(i)!= removeMask(i+1)) + rList(removeMask(i)) = i; + else + kList(i-removeMask(i)) = i; + } + ); + Kokkos::fence(); + + removeList = rList; + keepList = kList; + +} + + +void pFlow::boundaryBaseKernels::createRemoveKeepIndices +( + deviceViewType1D indices, + uint32 numRemove, + deviceViewType1D removeMask, + deviceViewType1D& removeIndices, + deviceViewType1D& keepIndices +) +{ + uint32 numTotal = indices.size(); + uint32 numKeep = numTotal - numRemove; + deviceViewType1D rIndices("rIndices",numRemove); + deviceViewType1D kIndices("kIndices",numKeep); + + exclusiveScan(removeMask, 0u, numTotal+1, removeMask, 0u); + + Kokkos::parallel_for + ( + "pFlow::boundaryBaseKernels::createRemoveKeepLists", + deviceRPolicyStatic(0, numTotal), + LAMBDA_HD(uint32 i) + { + if(removeMask(i)!= removeMask(i+1)) + rIndices(removeMask(i)) = indices(i); + else + kIndices(i-removeMask(i)) = indices(i); + } + ); + Kokkos::fence(); + + removeIndices = rIndices; + keepIndices = kIndices; +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.hpp new file mode 100644 index 00000000..e3166721 --- /dev/null +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBaseKernels.hpp @@ -0,0 +1,46 @@ +/*------------------------------- 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 __boundaryBaseKernels_hpp__ +#define __boundaryBaseKernels_hpp__ + +#include "phasicFlowKokkos.hpp" + +namespace pFlow::boundaryBaseKernels +{ + +void createRemoveKeepLists( + uint32 numTotal, + uint32 numRemove, + deviceViewType1D removeMask, + deviceViewType1D& removeList, + deviceViewType1D& keepList); + +void createRemoveKeepIndices( + deviceViewType1D indices, + uint32 numRemove, + deviceViewType1D removeMask, + deviceViewType1D& removeIndices, + deviceViewType1D& keepIndices); + + +} + +#endif //__boundaryBaseKernels_hpp__ \ No newline at end of file diff --git a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp index 9c0706d4..b48e80ca 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryExit/boundaryExit.cpp @@ -52,106 +52,36 @@ bool pFlow::boundaryExit::beforeIteration } uint32 s = size(); - deviceViewType1D delFlags("delFlags",s+1); - deviceViewType1D keepFlags("keepFlags", s+1); - fill(delFlags, 0, s+1, 0u); - fill(keepFlags, 0, s+1, 0u); - - using policy = Kokkos::RangePolicy< - pFlow::DefaultExecutionSpace, - Kokkos::Schedule, - Kokkos::IndexType>; + deviceViewType1D deleteFlags("deleteFlags",s+1); + fill(deleteFlags, 0, s+1, 0u); auto points = thisPoints(); - uint32 numDeleted = 0; auto p = boundaryPlane().infPlane(); - + + uint32 numDeleted = 0; Kokkos::parallel_reduce ( "boundaryExit::beforeIteration", - policy(0,s), + deviceRPolicyStatic(0,s), LAMBDA_HD(uint32 i, uint32& delToUpdate) { if(p.pointInNegativeSide(points(i))) { - delFlags(i)=1; + deleteFlags(i)=1; delToUpdate++; } - else - { - keepFlags(i) = 1; - } }, numDeleted ); - - + // no point is deleted if(numDeleted == 0 ) { return true; } - exclusiveScan(delFlags, 0u, s+1, delFlags, 0u); - exclusiveScan(keepFlags, 0u, s+1, keepFlags, 0u); - - deviceViewType1D deleteList("deleteList", numDeleted); - - Kokkos::parallel_for - ( - "boundaryExit::parllel_for", - policy(0, size()), - LAMBDA_HD(uint32 i) - { - if(delFlags(i)!= delFlags(i+1)) - deleteList(delFlags(i)) = i; - } - ); - Kokkos::fence(); - - deviceScatteredFieldAccess deleteIndices( - numDeleted, - deleteList, - indexList().deviceViewAll()); - - // tell internal to remove these points from its list - if(!internal().deletePoints(deleteIndices)) - { - fatalErrorInFunction<< - "error in deleting points from boundary "<< name()< newIndices("newIndices", newSize); - auto oldIndices = indexList().deviceViewAll(); - - Kokkos::parallel_for( - "fillIndices", - policy(0,s), - LAMBDA_HD(uint32 i) - { - if(keepFlags(i)!= keepFlags(i+1)) - { - newIndices(keepFlags(i)) = oldIndices(i); - } - } - ); - Kokkos::fence(); - setNewIndices(newIndices); - } - - //TODO: notify observers about changes - - return true; + return this->removeIndices(numDeleted, deleteFlags); } bool pFlow::boundaryExit::iterate diff --git a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp index 8cd8864f..2e6d20dc 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryList.cpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryList.cpp @@ -40,9 +40,19 @@ bool pFlow::boundaryList::updateLists() dist[4] = boundary(4).neighborLength(); dist[5] = boundary(5).neighborLength(); + realx3 lowerExt = + boundary(0).boundaryExtensionLength() + + boundary(2).boundaryExtensionLength() + + boundary(4).boundaryExtensionLength(); + realx3 upperExt = + boundary(1).boundaryExtensionLength()+ + boundary(3).boundaryExtensionLength()+ + boundary(5).boundaryExtensionLength(); + + auto extDomain = pStruct_.simDomain().extendThisDomain(lowerExt, upperExt); pStruct_.updateFlag( - pStruct_.simDomain().thisDomain(), + extDomain, dist); const auto& maskD = pStruct_.activePointsMaskDevice(); boundary(0).setSize( maskD.leftSize() ); diff --git a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp new file mode 100644 index 00000000..51fc7ab3 --- /dev/null +++ b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.cpp @@ -0,0 +1,125 @@ +/*------------------------------- 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 "boundaryPeriodic.hpp" +#include "internalPoints.hpp" +#include "dictionary.hpp" + + +pFlow::boundaryPeriodic::boundaryPeriodic +( + const dictionary& dict, + const plane& bplane, + internalPoints& internal +) +: + boundaryBase(dict, bplane, internal), + mirrorBoundaryIndex_(dict.getVal("mirrorBoundaryIndex")) +{ + extendedPlane_ = boundaryBase::boundaryPlane().parallelPlane(-boundaryBase::neighborLength()); +} + +pFlow::real pFlow::boundaryPeriodic::neighborLength() const +{ + return 2.0*boundaryBase::neighborLength(); +} + +pFlow::realx3 pFlow::boundaryPeriodic::boundaryExtensionLength() const +{ + return -neighborLength()*boundaryBase::boundaryPlane().normal(); +} + +const pFlow::plane &pFlow::boundaryPeriodic::boundaryPlane() const +{ + return extendedPlane_; +} + +bool pFlow::boundaryPeriodic::beforeIteration( + uint32 iterNum, + real t, + real dt) +{ + // nothing have to be done + if(empty()) + { + return true; + } + + uint32 s = size(); + deviceViewType1D transferFlags("transferFlags",s+1); + fill(transferFlags, 0, s+1, 0u); + + auto points = thisPoints(); + auto p = boundaryPlane().infPlane(); + uint32 numTransfered = 0; + + Kokkos::parallel_reduce + ( + "boundaryPeriodic::beforeIteration", + deviceRPolicyStatic(0u,s), + LAMBDA_HD(uint32 i, uint32& trnasToUpdate) + { + if(p.pointInNegativeSide(points(i))) + { + transferFlags(i)=1; + trnasToUpdate++; + } + }, + numTransfered + ); + + // no point to be transfered + if(numTransfered == 0 ) + { + return true; + } + + // to obtain the transfer vector + const auto& thisP = boundaryPlane(); + const auto& mirrorP = internal().boundary(mirrorBoundaryIndex_).boundaryPlane(); + realx3 transferVec = thisP.normal()*(thisP.d() + mirrorP.d()); + + return transferPoints + ( + numTransfered, + transferFlags, + mirrorBoundaryIndex_, + transferVec + ); + +} + +bool pFlow::boundaryPeriodic::iterate +( + uint32 iterNum, + real t +) +{ + return true; +} + +bool pFlow::boundaryPeriodic::afterIteration +( + uint32 iterNum, + real t +) +{ + return true; +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp new file mode 100644 index 00000000..01223885 --- /dev/null +++ b/src/phasicFlow/structuredData/boundaries/boundaryPeriodic/boundaryPeriodic.hpp @@ -0,0 +1,76 @@ +/*------------------------------- 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 __boundaryPeriodic_hpp__ +#define __boundaryPeriodic_hpp__ + + +#include "boundaryBase.hpp" + +namespace pFlow +{ + +class boundaryPeriodic +: + public boundaryBase +{ +private: + + uint32 mirrorBoundaryIndex_; + + plane extendedPlane_; + +public: + + TypeInfo("boundary"); + + boundaryPeriodic( + const dictionary& dict, + const plane& bplane, + internalPoints& internal); + + virtual + ~boundaryPeriodic() override= default; + + add_vCtor + ( + boundaryBase, + boundaryPeriodic, + dictionary + ); + + real neighborLength()const override; + + realx3 boundaryExtensionLength()const override; + + const plane& boundaryPlane()const override; + + bool beforeIteration(uint32 iterNum, real t, real dt) override; + + bool iterate(uint32 iterNum, real t) override; + + bool afterIteration(uint32 iterNum, real t) override; + + +}; + +} + +#endif \ No newline at end of file diff --git a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp index 5cc3c6fb..59ba7aa1 100644 --- a/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp +++ b/src/phasicFlow/structuredData/domain/regularSimulationDomain.cpp @@ -8,8 +8,8 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts() { auto& boundaries = this->subDict("boundaries"); - this->addDict("regularSimulationDomainBoundaries", boundaries); - auto& rbBoundaries = this->subDict("regularSimulationDomainBoundaries"); + this->addDict("regularBoundaries", boundaries); + auto& rbBoundaries = this->subDict("regularBoundaries"); real neighborLength = boundaries.getVal("neighborLength"); @@ -37,7 +37,25 @@ bool pFlow::regularSimulationDomain::createBoundaryDicts() " in dictionary "<< boundaries.globalName()<("type"); + uint32 mirrorIndex = 0; + if(bType == "periodic") + { + if(bName == bundaryName(0)) mirrorIndex = 1; + if(bName == bundaryName(1)) mirrorIndex = 0; + if(bName == bundaryName(2)) mirrorIndex = 3; + if(bName == bundaryName(3)) mirrorIndex = 2; + if(bName == bundaryName(4)) mirrorIndex = 5; + if(bName == bundaryName(5)) mirrorIndex = 4; + if(! bDict.addOrReplace("mirrorBoundaryIndex", mirrorIndex)) + { + fatalErrorInFunction<< + "canno add entry mirroBoundaryIndex to dictionary "<< + bDict.globalName()<subDict("regularSimulationDomainBoundaries"); + return this->subDict("regularBoundaries"); } bool pFlow::regularSimulationDomain::requiresDataTransfer()const diff --git a/src/phasicFlow/structuredData/domain/simulationDomain.cpp b/src/phasicFlow/structuredData/domain/simulationDomain.cpp index 232706c8..cf4201a6 100644 --- a/src/phasicFlow/structuredData/domain/simulationDomain.cpp +++ b/src/phasicFlow/structuredData/domain/simulationDomain.cpp @@ -41,8 +41,19 @@ pFlow::simulationDomain::simulationDomain(systemControl& control) } -pFlow::uniquePtr - pFlow::simulationDomain::create(systemControl& control) +pFlow::domain pFlow::simulationDomain::extendThisDomain +( + const realx3 &lowerPointExtension, + const realx3 &upperPointExtension +) const +{ + realx3 minP = thisDomain().domainBox().minPoint() + lowerPointExtension; + realx3 maxP = thisDomain().domainBox().maxPoint() + upperPointExtension; + return domain({minP, maxP}); +} + +pFlow::uniquePtr +pFlow::simulationDomain::create(systemControl &control) { word sType = angleBracketsNames( "simulationDomain", diff --git a/src/phasicFlow/structuredData/domain/simulationDomain.hpp b/src/phasicFlow/structuredData/domain/simulationDomain.hpp index 35718798..af49ee21 100644 --- a/src/phasicFlow/structuredData/domain/simulationDomain.hpp +++ b/src/phasicFlow/structuredData/domain/simulationDomain.hpp @@ -162,6 +162,11 @@ public: return thisDomain_; } + + domain extendThisDomain( + const realx3& lowerPointExtension, + const realx3& upperPointExtension)const; + inline const auto& globalBoundaryDict()const { diff --git a/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp b/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp index 3f727513..e9508076 100644 --- a/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp +++ b/src/phasicFlow/structuredData/pointStructure/internalPoints.cpp @@ -22,6 +22,7 @@ Licence: #include "internalPoints.hpp" #include "domain.hpp" #include "Vectors.hpp" +#include "internalPointsKernels.hpp" void pFlow::internalPoints::syncPFlag()const { @@ -32,6 +33,54 @@ void pFlow::internalPoints::syncPFlag()const } } +bool pFlow::internalPoints::deletePoints(deviceViewType1D delPoints) +{ + if(!pFlagsD_.deletePoints(delPoints)) + { + fatalErrorInFunction<< + "Error in deleting points from internal points"< changePoints, + uint32 boundaryIndex +) +{ + if(boundaryIndex>5) + { + fatalErrorInFunction<< + "Invalid boundary index "<< boundaryIndex< changePoints, + realx3 transferVector +) +{ + pFlow::internalPointsKernels::changePosition + ( + pointPosition_.deviceViewAll(), + changePoints, + transferVector + ); + WARNING<<"Notify about the change in the position of points"< delPoints); + + bool changePointsFlag( + deviceViewType1D changePoints, + uint32 boundaryIndex); + + bool changePointsPoisition( + deviceViewType1D changePoints, + realx3 transferVector); + public: //friend class dynamicinternalPoints; @@ -165,6 +180,18 @@ public: { return pFlagsD_.activeRange(); } + + virtual + Time& time() = 0; + + virtual + const Time& time() const = 0; + + virtual + boundaryBase& boundary(size_t boundaryIndex ) = 0; + + virtual + const boundaryBase& boundary(size_t boundaryIndex ) const = 0; ///@brief delete points at indices given in delPoints. /// The default is that delPoints contains sorted indices diff --git a/src/phasicFlow/structuredData/pointStructure/internalPointsKernels.cpp b/src/phasicFlow/structuredData/pointStructure/internalPointsKernels.cpp new file mode 100644 index 00000000..e5061d2a --- /dev/null +++ b/src/phasicFlow/structuredData/pointStructure/internalPointsKernels.cpp @@ -0,0 +1,44 @@ +/*------------------------------- 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 "internalPointsKernels.hpp" + +void pFlow::internalPointsKernels::changePosition +( + deviceViewType1D points, + deviceViewType1D indices, + realx3 transferVector +) +{ + auto s = static_cast(indices.size()); + if(s==0u)return; + + Kokkos::parallel_for + ( + "internalPointsKernels::changePosition", + deviceRPolicyStatic(0u, s), + LAMBDA_HD(uint32 i) + { + points[indices[i]] += transferVector; + } + ); + Kokkos::fence(); + +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/pointStructure/internalPointsKernels.hpp b/src/phasicFlow/structuredData/pointStructure/internalPointsKernels.hpp new file mode 100644 index 00000000..2b31f294 --- /dev/null +++ b/src/phasicFlow/structuredData/pointStructure/internalPointsKernels.hpp @@ -0,0 +1,31 @@ +/*------------------------------- 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 "phasicFlowKokkos.hpp" + +namespace pFlow::internalPointsKernels +{ + +void changePosition( + deviceViewType1D points, + deviceViewType1D indices, + realx3 transferVector); + +} \ No newline at end of file diff --git a/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp b/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp index f39b702d..ef91ccff 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp +++ b/src/phasicFlow/structuredData/pointStructure/pointFlag.hpp @@ -51,6 +51,10 @@ class pointFlag using device_type = typename viewType::device_type; + using rPolicy = Kokkos::RangePolicy< + execution_space, + Kokkos::IndexType>; + protected: viewType flags_; @@ -75,7 +79,34 @@ protected: //- Protected methods - + uint8 getBoundaryFlag(uint32 index)const + { + uint8 flg; + switch (index){ + case 0u: + flg = Flag::LEFT; + break; + case 1u: + flg = Flag::RIGHT; + break; + case 2u: + flg = Flag::BOTTOM; + break; + case 3u: + flg = Flag::TOP; + break; + case 4u: + flg = Flag::REAR; + break; + case 5u: + flg = Flag::FRONT; + break; + default: + flg=0; + } + + return flg; + } public: @@ -255,9 +286,15 @@ public: const box& validBox, ViewType1D points); - bool deletePoints( scatteredFieldAccess points); + + bool deletePoints( + ViewType1D points); + + bool changeFlags( + ViewType1D changePoints, + uint32 boundaryIndex); /// @brief mark points based on their position in the domain. /// This should be the first method to be called when updating diff --git a/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp b/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp index 88e3d145..aad7ff6e 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp +++ b/src/phasicFlow/structuredData/pointStructure/pointFlagKernels.hpp @@ -158,12 +158,59 @@ bool pFlow::pointFlag::deletePoints { if(points.empty())return true; - uint32 minIndex = points.getFirstCopy(); - uint32 maxIndex = points.getLastCopy(); + //uint32 minIndex = points.getFirstCopy(); + //uint32 maxIndex = points.getLastCopy(); - if(maxIndexactiveRange_.end())return false; - if(minIndexactiveRange_.end())return false; + //if(minIndex>; + + uint32 numDeleted = 0; + Kokkos::parallel_reduce + ( + "pointFlagKernels::deletePoints", + policy(0u, points.size()), + CLASS_LAMBDA_HD(uint32 i, uint32& valDelUpdate) + { + uint32 n = points(i); + if(isActive(n)) + { + valDelUpdate++; + flags_[n] = Flag::DELETED; + } + }, + numDeleted + ); + + if(numDeleted >= numActive_) + { + activeRange_ = {0, 0}; + numDeleted == numActive_; + } + + numActive_ = numActive_ - numDeleted; + isAllActive_ = + (activeRange_.numElements() == numActive_)&& numActive_>0; + + return true; +} +template +bool pFlow::pointFlag::deletePoints +( + ViewType1D points +) +{ + + uint32 s = points.size(); + if(s==0u)return true; + + //if(maxIndexactiveRange_.end())return false; + //if(minIndex::deletePoints } +template +bool pFlow::pointFlag::changeFlags +( + ViewType1D changePoints, + uint32 boundaryIndex +) +{ + auto flg = getBoundaryFlag(boundaryIndex); + Kokkos::parallel_for + ( + "pointFlag::changeFlags", + rPolicy(0, changePoints.size()), + LAMBDA_HD(uint32 i) + { + flags_[changePoints(i)] = flg; + } + ); + Kokkos::fence(); + return true; +} + template pFlow::uint32 pFlow::pointFlag::markPointRegions ( diff --git a/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp b/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp index 216d5585..be5b6998 100644 --- a/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp +++ b/src/phasicFlow/structuredData/pointStructure/pointStructure.hpp @@ -119,12 +119,23 @@ public: return boundaries_; } - auto& boundary(size_t i ) + + Time& time() override + { + return demComponent::time(); + } + + const Time& time() const override + { + return demComponent::time(); + } + + boundaryBase& boundary(size_t i ) override { return boundaries_[i]; } - auto& boundary(size_t i)const + const boundaryBase& boundary(size_t i)const override { return boundaries_[i]; } @@ -134,6 +145,7 @@ public: return simulationDomain_(); } + // - IO methods /// @brief read the point structure from the input stream.