diff --git a/CMakeLists.txt b/CMakeLists.txt index fa88bc95..9e770954 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -76,7 +76,7 @@ add_subdirectory(solvers) add_subdirectory(utilities) #add_subdirectory(DEMSystems) -add_subdirectory(testIO) + install(FILES "${PROJECT_BINARY_DIR}/phasicFlowConfig.H" DESTINATION include) diff --git a/cmake/bashrc b/cmake/bashrc index 32c92d87..d5e04cf8 100644 --- a/cmake/bashrc +++ b/cmake/bashrc @@ -1,7 +1,7 @@ export pFlow_PROJECT_VERSION=v-1.0 -export pFlow_PROJECT=phasicFlow +export pFlow_PROJECT=phasicFlowV-1.0 projectDir="$HOME/PhasicFlow" diff --git a/solvers/CMakeLists.txt b/solvers/CMakeLists.txt index 7db5fd36..5ecf3fb5 100644 --- a/solvers/CMakeLists.txt +++ b/solvers/CMakeLists.txt @@ -5,3 +5,5 @@ add_subdirectory(iterateSphereParticles) add_subdirectory(iterateGeometry) add_subdirectory(sphereGranFlow) + +add_subdirectory(grainGranFlow) diff --git a/solvers/grainGranFlow/CMakeLists.txt b/solvers/grainGranFlow/CMakeLists.txt new file mode 100755 index 00000000..afd35eb4 --- /dev/null +++ b/solvers/grainGranFlow/CMakeLists.txt @@ -0,0 +1,7 @@ + +set(source_files +grainGranFlow.cpp +) +set(link_lib Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities) + +pFlow_make_executable_install(grainGranFlow source_files link_lib) diff --git a/solvers/grainGranFlow/createDEMComponents.hpp b/solvers/grainGranFlow/createDEMComponents.hpp new file mode 100755 index 00000000..a92ed5c4 --- /dev/null +++ b/solvers/grainGranFlow/createDEMComponents.hpp @@ -0,0 +1,50 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ + +// +REPORT(0)<<"\nReading sphere particles . . ."<( + objectFile( + insertionFile__, + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS + ), + sphParticles, + sphParticles.shapes() + );*/ + +auto grnInsertion = pFlow::grainInsertion( + grnParticles, + grnParticles.grains()); + +REPORT(0)<<"\nCreating interaction model for sphere-sphere contact and sphere-wall contact . . ."< +class cGAbsoluteLinear +{ +public: + + struct contactForceStorage + { + realx3 overlap_t_ = 0.0; + }; + + + struct linearProperties + { + real kn_ = 1000.0; + real kt_ = 800.0; + real en_ = 0.0; + real ethat_ = 0.0; + real mu_ = 0.00001; + + INLINE_FUNCTION_HD + linearProperties(){} + + INLINE_FUNCTION_HD + linearProperties(real kn, real kt, real en, real etha_t, real mu ): + kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu) + {} + + INLINE_FUNCTION_HD + linearProperties(const linearProperties&)=default; + + INLINE_FUNCTION_HD + linearProperties& operator=(const linearProperties&)=default; + + INLINE_FUNCTION_HD + ~linearProperties() = default; + }; + + + +protected: + + using LinearArrayType = symArray; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + LinearArrayType linearProperties_; + + int32 addDissipationModel_; + + + bool readLinearDictionary(const dictionary& dict) + { + + + auto kn = dict.getVal("kn"); + auto kt = dict.getVal("kt"); + auto en = dict.getVal("en"); + auto et = dict.getVal("et"); + auto mu = dict.getVal("mu"); + + auto nElem = kn.size(); + + if(nElem != kt.size()) + { + fatalErrorInFunction<< + "sizes of kn("< prop("prop", nElem); + ForAll(i,kn) + { + prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] }; + } + + linearProperties_.assign(prop); + + auto adm = dict.getVal("additionalDissipationModel"); + + + if(adm == "none") + { + addDissipationModel_ = 1; + } + else if(adm == "LU") + { + addDissipationModel_ = 2; + } + else if (adm == "GB") + { + addDissipationModel_ = 3; + } + else + { + addDissipationModel_ = 1; + } + + return true; + + } + + static const char* modelName() + { + if constexpr (limited) + { + return "cGAbsoluteLinearLimited"; + } + else + { + return "cGAbsoluteLinearNonLimited"; + } + return ""; + } + +public: + + + TypeInfoNV(modelName()); + + INLINE_FUNCTION_HD + cGAbsoluteLinear(){} + + cGAbsoluteLinear(int32 nMaterial, const ViewType1D& rho, const dictionary& dict) + : + numMaterial_(nMaterial), + rho_("rho",nMaterial), + linearProperties_("linearProperties",nMaterial) + { + + Kokkos::deep_copy(rho_,rho); + if(!readLinearDictionary(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + cGAbsoluteLinear(const cGAbsoluteLinear&) = default; + + INLINE_FUNCTION_HD + cGAbsoluteLinear(cGAbsoluteLinear&&) = default; + + INLINE_FUNCTION_HD + cGAbsoluteLinear& operator=(const cGAbsoluteLinear&) = default; + + INLINE_FUNCTION_HD + cGAbsoluteLinear& operator=(cGAbsoluteLinear&&) = default; + + + INLINE_FUNCTION_HD + ~cGAbsoluteLinear()=default; + + INLINE_FUNCTION_HD + int32 numMaterial()const + { + return numMaterial_; + } + + //// - Methods + INLINE_FUNCTION_HD + void contactForce + ( + const real dt, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, + const real Ri, + const real Rj, + const real cGFi, + const real cGFj, + const real ovrlp_n, + const realx3& Vr, + const realx3& Nij, + contactForceStorage& history, + realx3& FCn, + realx3& FCt + )const + { + + + auto prop = linearProperties_(propId_i,propId_j); + + + real f_ = ( cGFi + cGFj )/2 ; + + + real vrn = dot(Vr, Nij); + realx3 Vt = Vr - vrn*Nij; + + history.overlap_t_ += Vt*dt; + + real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j]; + + real sqrt_meff = sqrt((mi*mj)/(mi+mj)); + + + // disipation model + if (addDissipationModel_==2) + { + prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_)); + } + else if (addDissipationModel_==3) + { + auto pie =3.14; + prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) )); + } + + + + + real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/ + sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0)); + + //REPORT(0)<<"\n en n is : "< ft_fric) + { + if( length(history.overlap_t_) >static_cast(0.0)) + { + if constexpr (limited) + { + FCt *= (ft_fric/ft); + history.overlap_t_ = - (FCt/prop.kt_); + } + else + { + FCt = (FCt/ft)*ft_fric; + } + //cout<<"friction is applied here \n"; + + } + else + { + FCt = 0.0; + } + } + + } + +}; + +} //pFlow::cfModels + +#endif diff --git a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp new file mode 100755 index 00000000..48a77884 --- /dev/null +++ b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp @@ -0,0 +1,340 @@ +/*------------------------------- 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 __cGRelativeLinearCF_hpp__ +#define __cGRelativeLinearCF_hpp__ + +#include "types.hpp" +#include "symArrays.hpp" + + + + + + + + +namespace pFlow::cfModels +{ + +template +class cGRelativeLinear +{ +public: + + struct contactForceStorage + { + realx3 overlap_t_ = 0.0; + }; + + + struct linearProperties + { + real kn_ = 1000.0; + real kt_ = 800.0; + real en_ = 0.0; + real ethat_ = 0.0; + real mu_ = 0.00001; + + INLINE_FUNCTION_HD + linearProperties(){} + + INLINE_FUNCTION_HD + linearProperties(real kn, real kt, real en, real etha_t, real mu ): + kn_(kn), kt_(kt), en_(en),ethat_(etha_t), mu_(mu) + {} + + INLINE_FUNCTION_HD + linearProperties(const linearProperties&)=default; + + INLINE_FUNCTION_HD + linearProperties& operator=(const linearProperties&)=default; + + INLINE_FUNCTION_HD + ~linearProperties() = default; + }; + + + +protected: + + using LinearArrayType = symArray; + + int32 numMaterial_ = 0; + + ViewType1D rho_; + + LinearArrayType linearProperties_; + + int32 addDissipationModel_; + + + bool readLinearDictionary(const dictionary& dict) + { + + + auto kn = dict.getVal("kn"); + auto kt = dict.getVal("kt"); + auto en = dict.getVal("en"); + auto et = dict.getVal("et"); + auto mu = dict.getVal("mu"); + + + auto nElem = kn.size(); + + + if(nElem != kt.size()) + { + fatalErrorInFunction<< + "sizes of kn("< prop("prop", nElem); + ForAll(i,kn) + { + prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] }; + } + + linearProperties_.assign(prop); + + auto adm = dict.getVal("additionalDissipationModel"); + + + if(adm == "none") + { + addDissipationModel_ = 1; + } + else if(adm == "LU") + { + addDissipationModel_ = 2; + } + else if (adm == "GB") + { + addDissipationModel_ = 3; + } + else + { + addDissipationModel_ = 1; + } + + return true; + + } + + static const char* modelName() + { + if constexpr (limited) + { + return "cGRelativeLinearLimited"; + } + else + { + return "cGRelativeLinearNonLimited"; + } + return ""; + } + +public: + + + TypeInfoNV(modelName()); + + INLINE_FUNCTION_HD + cGRelativeLinear(){} + + cGRelativeLinear(int32 nMaterial, const ViewType1D& rho, const dictionary& dict) + : + numMaterial_(nMaterial), + rho_("rho",nMaterial), + linearProperties_("linearProperties",nMaterial) + { + + Kokkos::deep_copy(rho_,rho); + if(!readLinearDictionary(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + cGRelativeLinear(const cGRelativeLinear&) = default; + + INLINE_FUNCTION_HD + cGRelativeLinear(cGRelativeLinear&&) = default; + + INLINE_FUNCTION_HD + cGRelativeLinear& operator=(const cGRelativeLinear&) = default; + + INLINE_FUNCTION_HD + cGRelativeLinear& operator=(cGRelativeLinear&&) = default; + + + INLINE_FUNCTION_HD + ~cGRelativeLinear()=default; + + INLINE_FUNCTION_HD + int32 numMaterial()const + { + return numMaterial_; + } + + //// - Methods + INLINE_FUNCTION_HD + void contactForce + ( + const real dt, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, + const real Ri, + const real Rj, + const real cGFi, + const real cGFj, + const real ovrlp_n, + const realx3& Vr, + const realx3& Nij, + contactForceStorage& history, + realx3& FCn, + realx3& FCt + )const + { + + auto prop = linearProperties_(propId_i,propId_j); + + real f_ = ( cGFi + cGFj )/2 ; + + + real vrn = dot(Vr, Nij); + realx3 Vt = Vr - vrn*Nij; + + history.overlap_t_ += Vt*dt; + + real mi = 3*Pi/4*pow(Ri,3.0)*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,3.0)*rho_[propId_j]; + + real sqrt_meff = sqrt((mi*mj)/(mi+mj)); + + + if (addDissipationModel_==2) + { + prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_)); + } + else if (addDissipationModel_==3) + { + auto pie =3.14; + prop.en_ = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(pie,2)))) ) )); + } + real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/ + sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0)); + + + FCn = ( -pow(f_,1.0)*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; + FCt = ( -pow(f_,1.0)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,0.5) * prop.ethat_*Vt); + + + + real ft = length(FCt); + real ft_fric = prop.mu_ * length(FCn); + + if(ft > ft_fric) + { + if( length(history.overlap_t_) >static_cast(0.0)) + { + if constexpr (limited) + { + FCt *= (ft_fric/ft); + history.overlap_t_ = - (FCt/prop.kt_); + } + else + { + FCt = (FCt/ft)*ft_fric; + } + //cout<<"friction is applied here \n"; + + } + else + { + FCt = 0.0; + } + } + + } + +}; + +} //pFlow::cfModels + +#endif diff --git a/src/Interaction/Models/grainContactForceModels.hpp b/src/Interaction/Models/grainContactForceModels.hpp new file mode 100755 index 00000000..908117f5 --- /dev/null +++ b/src/Interaction/Models/grainContactForceModels.hpp @@ -0,0 +1,45 @@ +/*------------------------------- 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 __grainContactForceModels_hpp__ +#define __grainContactForceModels_hpp__ + +#include "cGAbsoluteLinearCF.hpp" +#include "cGRelativeLinearCF.hpp" + +#include "grainRolling.hpp" + + + +namespace pFlow::cfModels +{ + + +using limitedCGAbsoluteLinearGrainRolling = grainRolling>; +using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling>; + +using limitedCGRelativeLinearGrainRolling = grainRolling>; +using nonLimitedCGRelativeLinearGrainRolling = grainRolling>; + +} + + + +#endif //__grainContactForceModels_hpp__ diff --git a/src/Interaction/Models/rolling/grainRolling.hpp b/src/Interaction/Models/rolling/grainRolling.hpp new file mode 100755 index 00000000..b8eff130 --- /dev/null +++ b/src/Interaction/Models/rolling/grainRolling.hpp @@ -0,0 +1,119 @@ +/*------------------------------- 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 __grainRolling_hpp__ +#define __grainRolling_hpp__ + + +namespace pFlow::cfModels +{ + +template +class grainRolling +: + public contactForceModel +{ +public: + + using contactForceStorage = + typename contactForceModel::contactForceStorage; + + + realSymArray_D mur_; + + bool readGrainDict(const dictionary& dict) + { + auto mur = dict.getVal("mur"); + + uint32 nMat; + + if(!realSymArray_D::getN(mur.size(),nMat) || nMat != this->numMaterial()) + { + fatalErrorInFunction<< + "wrong number of values supplied in mur. \n"; + return false; + } + + mur_.assign(mur); + + return true; + } + +public: + + TypeInfoNV(word("normal<"+contactForceModel::TYPENAME()+">")); + + + grainRolling(int32 nMaterial, const ViewType1D& rho, const dictionary& dict) + : + contactForceModel(nMaterial, rho, dict), + mur_("mur", nMaterial) + { + if(!readGrainDict(dict)) + { + fatalExit; + } + } + + INLINE_FUNCTION_HD + void rollingFriction + ( + const real dt, + const uint32 i, + const uint32 j, + const uint32 propId_i, + const uint32 propId_j, + const real Ri, + const real Rj, + const real cGFi, + const real cGFj, + const realx3& wi, + const realx3& wj, + const realx3& Nij, + const realx3& FCn, + realx3& Mri, + realx3& Mrj + )const + { + + realx3 w_hat = wi-wj; + real w_hat_mag = length(w_hat); + + if( !equal(w_hat_mag,0.0) ) + w_hat /= w_hat_mag; + else + w_hat = 0.0; + + auto Reff = (Ri*Rj)/(Ri+Rj); + + Mri = ( -mur_(propId_i,propId_j) *length(FCn) * Reff ) * w_hat ; + + //removing the normal part + // Mri = Mri - ( (Mri .dot. nij)*nij ) + + Mrj = -Mri; + } + + +}; + +} + +#endif diff --git a/src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.cpp b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.cpp new file mode 100644 index 00000000..a2949227 --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.cpp @@ -0,0 +1,104 @@ +#include "boundarySphereInteraction.hpp" +/*------------------------------- 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 +void pFlow::boundaryGrainInteraction::allocatePPPairs() +{ + ppPairs_.reset(nullptr); + ppPairs_ = makeUnique(1); +} + +template +void pFlow::boundaryGrainInteraction::allocatePWPairs() +{ + pwPairs_.reset(nullptr); + pwPairs_ = makeUnique(1); +} + + +template +pFlow::boundaryGrainInteraction::boundaryGrainInteraction( + const boundaryBase &boundary, + const grainParticles &grnPrtcls, + const GeometryMotionModel &geomMotion) + : generalBoundary(boundary, grnPrtcls.pStruct(), "", ""), + geometryMotion_(geomMotion), + grnParticles_(grnPrtcls) +{ +} + +template +pFlow::uniquePtr> +pFlow::boundaryGrainInteraction::create( + const boundaryBase &boundary, + const grainParticles &grnPrtcls, + const GeometryMotionModel &geomMotion) +{ + word cfTypeName = ContactForceModel::TYPENAME(); + word gmTypeName = MotionModel::TYPENAME(); + word bType = boundary.type(); + + word boundaryTypeName = angleBracketsNames3( + "boundaryGrainInteraction", + bType, + cfTypeName, + gmTypeName); + + word altBTypeName = angleBracketsNames2( + "boundaryGrainInteraction", + cfTypeName, + gmTypeName); + + if (boundaryBasevCtorSelector_.search(boundaryTypeName)) + { + pOutput.space(4) << "Creating boundry type " << Green_Text(boundaryTypeName) << + " for boundary " << boundary.name() << " . . ." << END_REPORT; + return boundaryBasevCtorSelector_[boundaryTypeName]( + boundary, + grnPrtcls, + geomMotion); + } + else if(boundaryBasevCtorSelector_[altBTypeName]) + { + // if boundary condition is not implemented, the default is used + + pOutput.space(4) << "Creating boundry type " << Green_Text(altBTypeName) << + " for boundary " << boundary.name() << " . . ." << END_REPORT; + return boundaryBasevCtorSelector_[altBTypeName]( + boundary, + grnPrtcls, + geomMotion); + } + else + { + printKeys + ( + fatalError << "Ctor Selector "<< boundaryTypeName<< + " and "<< altBTypeName << " do not exist. \n" + <<"Avaiable ones are: \n\n" + , + boundaryBasevCtorSelector_ + ); + fatalExit; + } + + return nullptr; +} diff --git a/src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.hpp b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.hpp new file mode 100644 index 00000000..4800194c --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteraction/boundaryGrainInteraction.hpp @@ -0,0 +1,186 @@ +/*------------------------------- 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 __boundaryGrainInteraction_hpp__ +#define __boundaryGrainInteraction_hpp__ + +#include "virtualConstructor.hpp" +#include "generalBoundary.hpp" +#include "sortedContactList.hpp" +#include "grainParticles.hpp" + +namespace pFlow +{ + +template +class boundaryGrainInteraction +: + public generalBoundary +{ +public: + + using BoundaryGrainInteractionType + = boundaryGrainInteraction; + + using GeometryMotionModel = geometryMotionModel; + + using ContactForceModel = contactForceModel; + + using MotionModel = typename geometryMotionModel::MotionModel; + + using ModelStorage = typename ContactForceModel::contactForceStorage; + + using IdType = uint32; + + using IndexType = uint32; + + using ContactListType = + sortedContactList; + +private: + + const GeometryMotionModel& geometryMotion_; + + /// const reference to sphere particles + const grainParticles& grnParticles_; + + uniquePtr ppPairs_ = nullptr; + + uniquePtr pwPairs_ = nullptr; + +protected: + + void allocatePPPairs(); + + void allocatePWPairs(); + +public: + + TypeInfoTemplate12("boundaryGrainInteraction", ContactForceModel, MotionModel); + + boundaryGrainInteraction( + const boundaryBase& boundary, + const grainParticles& grnPrtcls, + const GeometryMotionModel& geomMotion); + + create_vCtor + ( + BoundaryGrainInteractionType, + boundaryBase, + ( + const boundaryBase& boundary, + const grainParticles& grnPrtcls, + const GeometryMotionModel& geomMotion + ), + (boundary, grnPrtcls, geomMotion) + ); + + add_vCtor + ( + BoundaryGrainInteractionType, + BoundaryGrainInteractionType, + boundaryBase + ); + + ~boundaryGrainInteraction()override=default; + + const auto& grnParticles()const + { + return grnParticles_; + } + + const auto& geometryMotion()const + { + return geometryMotion_; + } + + ContactListType& ppPairs() + { + return ppPairs_(); + } + + const ContactListType& ppPairs()const + { + return ppPairs_(); + } + + ContactListType& pwPairs() + { + return pwPairs_(); + } + + const ContactListType& pwPairs()const + { + return pwPairs_(); + } + + bool ppPairsAllocated()const + { + if( ppPairs_)return true; + return false; + } + + bool pwPairsAllocated()const + { + if( pwPairs_)return true; + return false; + } + + virtual + bool grainGrainInteraction( + real dt, + const ContactForceModel& cfModel, + uint32 step) + { + // for default boundary, no thing to be done + return false; + } + + + + bool hearChanges + ( + real t, + real dt, + uint32 iter, + const message& msg, + const anyList& varList + ) override + { + + pOutput<<"Function (hearChanges in boundarySphereInteractions)is not implmented Message "<< + msg <boundaryName() <<" type "<< this->type()< create( + const boundaryBase& boundary, + const grainParticles& sphPrtcls, + const GeometryMotionModel& geomMotion + ); + +}; + +} + +#include "boundaryGrainInteraction.cpp" + +#endif //__boundaryGrainInteraction_hpp__ diff --git a/src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.cpp b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.cpp new file mode 100644 index 00000000..7875210d --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.cpp @@ -0,0 +1,23 @@ + + +template +pFlow::boundaryGrainInteractionList::boundaryGrainInteractionList +( + const grainParticles &grnPrtcls, + const gMModel &geomMotion +) +: + ListPtr>(6), + boundaries_(grnPrtcls.pStruct().boundaries()) +{ + //gSettings::sleepMiliSeconds(1000*pFlowProcessors().localRank()); + for(uint32 i=0; i<6; i++) + { + this->set( + i, + boundaryGrainInteraction::create( + boundaries_[i], + grnPrtcls, + geomMotion)); + } +} \ No newline at end of file diff --git a/src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.hpp b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.hpp new file mode 100644 index 00000000..ed4cabc8 --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/boundaryGrainInteractionList.hpp @@ -0,0 +1,40 @@ +#ifndef __boundaryGrainInteractionList_hpp__ +#define __boundaryGrainInteractionList_hpp__ + + +#include "boundaryList.hpp" +#include "ListPtr.hpp" +#include "boundaryGrainInteraction.hpp" + + +namespace pFlow +{ + + +template +class boundaryGrainInteractionList +: + public ListPtr> +{ +private: + + const boundaryList& boundaries_; + +public: + + boundaryGrainInteractionList( + const grainParticles& grnPrtcls, + const geometryMotionModel& geomMotion + ); + + ~boundaryGrainInteractionList()=default; + +}; + + + +} + +#include "boundaryGrainInteractionList.cpp" + +#endif //__boundaryGrainInteractionList_hpp__ \ No newline at end of file diff --git a/src/Interaction/grainInteraction/boundaries/createBoundaryGrainInteraction.hpp b/src/Interaction/grainInteraction/boundaries/createBoundaryGrainInteraction.hpp new file mode 100644 index 00000000..12049af0 --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/createBoundaryGrainInteraction.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 "boundaryGrainInteraction.hpp" +#include "periodicBoundaryGrainInteraction.hpp" + +#define createBoundaryGrainInteraction(ForceModel,GeomModel) \ + template class pFlow::boundaryGrainInteraction< \ + ForceModel, \ + GeomModel>; \ + \ + template class pFlow::periodicBoundaryGrainInteraction< \ + ForceModel, \ + GeomModel>; + diff --git a/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.cpp b/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.cpp new file mode 100644 index 00000000..0eedf973 --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.cpp @@ -0,0 +1,70 @@ +/*------------------------------- 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 "periodicBoundarySIKernels.hpp" + +template +pFlow::periodicBoundaryGrainInteraction::periodicBoundaryGrainInteraction( + const boundaryBase &boundary, + const grainParticles &grnPrtcls, + const GeometryMotionModel &geomMotion) + : boundaryGrainInteraction(boundary, grnPrtcls, geomMotion), + transferVec_(boundary.mirrorBoundary().displacementVectroToMirror()) +{ + if(boundary.thisBoundaryIndex()%2==1) + { + masterInteraction_ = true; + this->allocatePPPairs(); + this->allocatePWPairs(); + + } + else + { + masterInteraction_ = false; + } +} + +template +bool pFlow::periodicBoundaryGrainInteraction::grainGrainInteraction +( + real dt, + const ContactForceModel &cfModel, + uint32 step +) +{ + if(!masterInteraction_) return false; + + pFlow::periodicBoundarySIKernels::grainGrainInteraction( + dt, + this->ppPairs(), + cfModel, + transferVec_, + this->boundary().thisPoints(), + this->mirrorBoundary().thisPoints(), + this->grnParticles().diameter().deviceViewAll(), + this->grnParticles().coarseGrainFactor().deviceViewAll(), + this->grnParticles().propertyId().deviceViewAll(), + this->grnParticles().velocity().deviceViewAll(), + this->grnParticles().rVelocity().deviceViewAll(), + this->grnParticles().contactForce().deviceViewAll(), + this->grnParticles().contactTorque().deviceViewAll()); + + return false; +} \ No newline at end of file diff --git a/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.hpp b/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.hpp new file mode 100644 index 00000000..e9148c9d --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundaryGrainInteraction.hpp @@ -0,0 +1,96 @@ +/*------------------------------- 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 __periodicBoundaryGrainInteraction_hpp__ +#define __periodicBoundaryGrainInteraction_hpp__ + +#include "boundaryGrainInteraction.hpp" + +namespace pFlow +{ + +template +class periodicBoundaryGrainInteraction +: + public boundaryGrainInteraction +{ +public: + + using PBSInteractionType = + periodicBoundaryGrainInteraction; + + using BSInteractionType = + boundaryGrainInteraction; + + 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: + + realx3 transferVec_; + + bool masterInteraction_; +public: + + TypeInfoTemplate22("boundaryGrainInteraction", "periodic",ContactForceModel, MotionModel); + + + periodicBoundaryGrainInteraction( + const boundaryBase& boundary, + const grainParticles& grnPrtcls, + const GeometryMotionModel& geomMotion + ); + + add_vCtor + ( + BSInteractionType, + PBSInteractionType, + boundaryBase + ); + + ~periodicBoundaryGrainInteraction()override = default; + + + + + bool grainGrainInteraction( + real dt, + const ContactForceModel& cfModel, + uint32 step)override; + +}; + +} + +#include "periodicBoundaryGrainInteraction.cpp" + + +#endif //__periodicBoundaryGrainInteraction_hpp__ \ No newline at end of file diff --git a/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundarySIKernels.hpp b/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundarySIKernels.hpp new file mode 100644 index 00000000..6a8a618e --- /dev/null +++ b/src/Interaction/grainInteraction/boundaries/periodicBoundaryGrainInteraction/periodicBoundarySIKernels.hpp @@ -0,0 +1,121 @@ + +namespace pFlow::periodicBoundarySIKernels +{ + +template +inline +void grainGrainInteraction +( + real dt, + const ContactListType& cntctList, + const ContactForceModel& forceModel, + const realx3& transferVec, + const deviceScatteredFieldAccess& thisPoints, + const deviceScatteredFieldAccess& mirrorPoints, + const deviceViewType1D& diam, + const deviceViewType1D& coarseGrainFactor, + const deviceViewType1D& propId, + const deviceViewType1D& vel, + const deviceViewType1D& rVel, + const deviceViewType1D& cForce, + const deviceViewType1D& cTorque +) +{ + + using ValueType = typename ContactListType::ValueType; + uint32 ss = cntctList.size(); + uint32 lastItem = cntctList.loopCount(); + if(lastItem == 0u)return; + + Kokkos::parallel_for( + "pFlow::periodicBoundarySIKernels::grainGrainInteraction", + 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 = mirrorPoints.index(j); + + real Ri = 0.5*diam[ind_i]; + real Rj = 0.5*diam[ind_j]; + real cGFi = coarseGrainFactor[ind_i]; + real cGFj = coarseGrainFactor[ind_j]; + realx3 xi = thisPoints.field()[ind_i]; + realx3 xj = mirrorPoints.field()[ind_j]+transferVec; + real dist = length(xj-xi); + real ovrlp = (Ri+Rj) - dist; + + if( ovrlp >0.0 ) + { + auto Nij = (xj-xi)/dist; + auto wi = rVel[ind_i]; + auto wj = rVel[ind_j]; + auto Vr = vel[ind_i] - vel[ind_j] + cross((Ri*wi+Rj*wj), Nij); + + auto history = cntctList.getValue(n); + + int32 propId_i = propId[ind_i]; + int32 propId_j = propId[ind_j]; + + realx3 FCn, FCt, Mri, Mrj, Mij, Mji; + + // calculates contact force + forceModel.contactForce( + dt, i, j, + propId_i, propId_j, + Ri, Rj, cGFi , cGFj , + ovrlp, + Vr, Nij, + history, + FCn, FCt); + + forceModel.rollingFriction( + dt, i, j, + propId_i, propId_j, + Ri, Rj, cGFi , cGFj , + 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(&cForce[ind_i].x_,FC.x_); + Kokkos::atomic_add(&cForce[ind_i].y_,FC.y_); + Kokkos::atomic_add(&cForce[ind_i].z_,FC.z_); + + Kokkos::atomic_add(&cForce[ind_j].x_,-FC.x_); + Kokkos::atomic_add(&cForce[ind_j].y_,-FC.y_); + Kokkos::atomic_add(&cForce[ind_j].z_,-FC.z_); + + Kokkos::atomic_add(&cTorque[ind_i].x_, Mij.x_); + Kokkos::atomic_add(&cTorque[ind_i].y_, Mij.y_); + Kokkos::atomic_add(&cTorque[ind_i].z_, Mij.z_); + + Kokkos::atomic_add(&cTorque[ind_j].x_, Mji.x_); + Kokkos::atomic_add(&cTorque[ind_j].y_, Mji.y_); + Kokkos::atomic_add(&cTorque[ind_j].z_, Mji.z_); + + + cntctList.setValue(n,history); + + } + else + { + cntctList.setValue(n, ValueType()); + } + + }); + Kokkos::fence(); +} + + +} //pFlow::periodicBoundarySIKernels \ No newline at end of file diff --git a/src/Interaction/grainInteraction/grainInteraction/grainInteraction.cpp b/src/Interaction/grainInteraction/grainInteraction/grainInteraction.cpp new file mode 100644 index 00000000..65c84a77 --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteraction/grainInteraction.cpp @@ -0,0 +1,359 @@ +/*------------------------------- 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 class cLT> +bool pFlow::grainInteraction::createGrainInteraction() +{ + + realVector_D rhoD("densities", this->densities()); + + auto modelDict = this->subDict("model"); + + forceModel_ = makeUnique( + this->numMaterials(), + rhoD.deviceView(), + modelDict ); + + + uint32 nPrtcl = grnParticles_.size(); + + contactSearch_ = contactSearch::create( + subDict("contactSearch"), + grnParticles_.extendedDomain().domainBox(), + grnParticles_, + geometryMotion_, + timers()); + + ppContactList_ = makeUnique(nPrtcl+1); + + pwContactList_ = makeUnique(nPrtcl/5+1); + + return true; +} + + +template class cLT> +bool pFlow::grainInteraction::grainGrainInteraction() +{ + auto lastItem = ppContactList_().loopCount(); + + // create the kernel functor + pFlow::grainInteractionKernels::ppInteractionFunctor + ppInteraction( + this->dt(), + this->forceModel_(), + ppContactList_(), // to be read + grnParticles_.diameter().deviceViewAll(), + grnParticles_.coarseGrainFactor().deviceViewAll(), + grnParticles_.propertyId().deviceViewAll(), + grnParticles_.pointPosition().deviceViewAll(), + grnParticles_.velocity().deviceViewAll(), + grnParticles_.rVelocity().deviceViewAll(), + grnParticles_.contactForce().deviceViewAll(), + grnParticles_.contactTorque().deviceViewAll() + ); + + Kokkos::parallel_for( + "ppInteraction", + rpPPInteraction(0,lastItem), + ppInteraction + ); + + Kokkos::fence(); + + return true; +} + + +template class cLT> +bool pFlow::grainInteraction::grainWallInteraction() +{ + + uint32 lastItem = pwContactList_().loopCount(); + uint32 iter = this->currentIter(); + real t = this->currentTime(); + real dt = this->dt(); + + pFlow::grainInteractionKernels::pwInteractionFunctor + pwInteraction( + dt, + this->forceModel_(), + pwContactList_(), + geometryMotion_.getTriangleAccessor(), + geometryMotion_.getModel(iter, t, dt) , + grnParticles_.diameter().deviceViewAll() , + grnParticles_.coarseGrainFactor().deviceViewAll() , + grnParticles_.propertyId().deviceViewAll(), + grnParticles_.pointPosition().deviceViewAll(), + grnParticles_.velocity().deviceViewAll(), + grnParticles_.rVelocity().deviceViewAll() , + grnParticles_.contactForce().deviceViewAll(), + grnParticles_.contactTorque().deviceViewAll() , + geometryMotion_.triMotionIndex().deviceViewAll(), + geometryMotion_.propertyId().deviceViewAll(), + geometryMotion_.contactForceWall().deviceViewAll() + ); + + Kokkos::parallel_for( + "", + rpPWInteraction(0,lastItem), + pwInteraction + ); + + + Kokkos::fence(); + + return true; +} + + +template class cLT> +pFlow::grainInteraction::grainInteraction +( + systemControl& control, + const particles& prtcl, + const geometry& geom +) +: + interaction(control, prtcl, geom), + geometryMotion_(dynamic_cast(geom)), + grnParticles_(dynamic_cast(prtcl)), + boundaryInteraction_(grnParticles_,geometryMotion_), + ppInteractionTimer_("grain-graine interaction (internal)", &this->timers()), + pwInteractionTimer_("grain-wall interaction (internal)", &this->timers()), + boundaryInteractionTimer_("grain-grain interaction (boundary)",&this->timers()), + contactListMangementTimer_("list management (interal)", &this->timers()), + contactListMangementBoundaryTimer_("list management (boundary)", &this->timers()) +{ + + if(!createGrainInteraction()) + { + fatalExit; + } + + for(uint32 i=0; i<6; i++) + { + activeBoundaries_[i] = boundaryInteraction_[i].ppPairsAllocated(); + } +} + +template class cLT> +bool pFlow::grainInteraction::beforeIteration() +{ + return true; +} + +template class cLT> +bool pFlow::grainInteraction::iterate() +{ + + timeInfo ti = this->TimeInfo(); + auto iter = ti.iter(); + auto t = ti.t(); + auto dt = ti.dt(); + + auto& contactSearchRef = contactSearch_(); + + bool broadSearch = contactSearchRef.enterBroadSearch(ti); + bool broadSearchBoundary = contactSearchRef.enterBroadSearchBoundary(ti); + + /// update boundaries of the fields that are being used by inreaction + grnParticles_.diameter().updateBoundaries(DataDirection::SlaveToMaster); + grnParticles_.velocity().updateBoundaries(DataDirection::SlaveToMaster); + grnParticles_.rVelocity().updateBoundaries(DataDirection::SlaveToMaster); + grnParticles_.propertyId().updateBoundaries(DataDirection::SlaveToMaster); + + /// lists + if(broadSearch) + { + contactListMangementTimer_.start(); + ComputationTimer().start(); + ppContactList_().beforeBroadSearch(); + pwContactList_().beforeBroadSearch(); + ComputationTimer().end(); + contactListMangementTimer_.pause(); + } + + if(broadSearchBoundary) + { + contactListMangementBoundaryTimer_.start(); + ComputationTimer().start(); + for(uint32 i=0; i<6u; i++) + { + if(activeBoundaries_[i]) + { + auto& BI = boundaryInteraction_[i]; + BI.ppPairs().beforeBroadSearch(); + BI.pwPairs().beforeBroadSearch(); + } + } + ComputationTimer().end(); + contactListMangementBoundaryTimer_.pause(); + } + + if( grnParticles_.numActive()<=0)return true; + + ComputationTimer().start(); + if( !contactSearchRef.broadSearch( + ti, + ppContactList_(), + pwContactList_()) ) + { + fatalErrorInFunction<< + "unable to perform broadSearch.\n"; + fatalExit; + } + + for(uint32 i=0; i<6u; i++) + { + if(activeBoundaries_[i]) + { + auto& BI = boundaryInteraction_[i]; + if(!contactSearchRef.boundaryBroadSearch( + i, + ti, + BI.ppPairs(), + BI.pwPairs()) + ) + { + fatalErrorInFunction<< + "failed to perform broadSearch for boundary index "<(true); + const auto& cfModel = this->forceModel_(); + uint32 step=1; + boundaryInteractionTimer_.start(); + ComputationTimer().start(); + while(requireStep.anyElement(true) && step <= 10) + { + for(uint32 i=0; i<6u; i++) + { + if(requireStep[i] ) + { + requireStep[i] = boundaryInteraction_[i].grainGrainInteraction( + dt, + this->forceModel_(), + step + ); + } + } + step++; + } + ComputationTimer().end(); + boundaryInteractionTimer_.pause(); + + + ppInteractionTimer_.start(); + ComputationTimer().start(); + grainGrainInteraction(); + ComputationTimer().end(); + ppInteractionTimer_.end(); + + + pwInteractionTimer_.start(); + ComputationTimer().start(); + grainWallInteraction(); + ComputationTimer().start(); + pwInteractionTimer_.end(); + + { + boundaryInteractionTimer_.resume(); + ComputationTimer().start(); + auto requireStep = boundariesMask<6>(true); + + uint32 step = 11; + const auto& cfModel = this->forceModel_(); + while( requireStep.anyElement(true) && step < 20 ) + { + for(uint32 i=0; i<6u; i++) + { + if(requireStep[i]) + { + requireStep[i] = boundaryInteraction_[i].grainGrainInteraction( + dt, + cfModel, + step + ); + } + } + step++; + } + ComputationTimer().end(); + boundaryInteractionTimer_.end(); + } + + return true; +} + +template class cLT> +bool pFlow::grainInteraction::afterIteration() +{ + return true; +} + +template class cLT> +bool pFlow::grainInteraction::hearChanges +( + real t, + real dt, + uint32 iter, + const message& msg, + const anyList& varList +) +{ + if(msg.equivalentTo(message::ITEM_REARRANGE)) + { + notImplementedFunction; + } + return true; +} \ No newline at end of file diff --git a/src/Interaction/grainInteraction/grainInteraction/grainInteraction.hpp b/src/Interaction/grainInteraction/grainInteraction/grainInteraction.hpp new file mode 100644 index 00000000..2d30165b --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteraction/grainInteraction.hpp @@ -0,0 +1,169 @@ +/*------------------------------- 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 __grainInteraction_hpp__ +#define __grainInteraction_hpp__ + +#include "interaction.hpp" +#include "grainParticles.hpp" +#include "boundaryGrainInteractionList.hpp" +#include "grainInteractionKernels.hpp" +#include "boundariesMask.hpp" +#include "MPITimer.hpp" +//#include "unsortedContactList.hpp" + + + +namespace pFlow +{ + +template< + typename contactForceModel, + typename geometryMotionModel, + template class contactListType> +class grainInteraction +: + public interaction +{ +public: + + using GeometryMotionModel = geometryMotionModel; + + using ContactForceModel = contactForceModel; + + using MotionModel = typename geometryMotionModel::MotionModel; + + using ModelStorage = typename ContactForceModel::contactForceStorage; + + using BoundaryListType = boundaryGrainInteractionList; + + using IdType = uint32; + + using IndexType = uint32; + + using ContactListType = + contactListType; + + //using BoundaryContactListType = unsortedContactList; + + + +private: + + /// const reference to geometry + const GeometryMotionModel& geometryMotion_; + + /// const reference to particles + const grainParticles& grnParticles_; + + /// particle-particle and particle-wall interactions at boundaries + BoundaryListType boundaryInteraction_; + + /// a mask for active boundaries (boundaries with intreaction) + boundariesMask<6> activeBoundaries_; + + /// contact search object for pp and pw interactions + uniquePtr contactSearch_ = nullptr; + + /// contact force model + uniquePtr forceModel_ = nullptr; + + /// contact list for particle-particle interactoins (keeps the history) + uniquePtr ppContactList_ = nullptr; + + /// contact list for particle-wall interactions (keeps the history) + uniquePtr pwContactList_ = nullptr; + + + /// timer for particle-particle interaction computations + Timer ppInteractionTimer_; + + /// timer for particle-wall interaction computations + Timer pwInteractionTimer_; + + /// timer for boundary interaction time + Timer boundaryInteractionTimer_; + + /// timer for managing contact lists (only inernal points) + Timer contactListMangementTimer_; + + Timer contactListMangementBoundaryTimer_; + + + + bool createGrainInteraction(); + + bool grainGrainInteraction(); + + bool grainWallInteraction(); + + /// range policy for p-p interaction execution + using rpPPInteraction = + Kokkos::RangePolicy, Kokkos::Schedule>; + + /// range policy for p-w interaction execution + using rpPWInteraction = rpPPInteraction; + +public: + + TypeInfoTemplate13("grainInteraction", ContactForceModel, MotionModel, ContactListType); + + /// Constructor from components + grainInteraction( + systemControl& control, + const particles& prtcl, + const geometry& geom); + + + /// Add virtual constructor + add_vCtor + ( + interaction, + grainInteraction, + systemControl + ); + + /// This is called in time loop, before iterate. (overriden from demComponent) + bool beforeIteration() override; + + /// This is called in time loop. Perform the main calculations + /// when the component should evolve along time. (overriden from demComponent) + bool iterate() override; + + /// This is called in time loop, after iterate. (overriden from demComponent) + bool afterIteration() override; + + /// Check for changes in the point structures. (overriden from observer) + bool hearChanges( + real t, + real dt, + uint32 iter, + const message& msg, + const anyList& varList)override; + + +}; + + +} + +#include "grainInteraction.cpp" + +#endif //__grainInteraction_hpp__ diff --git a/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp b/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp new file mode 100644 index 00000000..7ee5ef39 --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteraction/grainInteractionKernels.hpp @@ -0,0 +1,337 @@ +/*------------------------------- 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 __grainInteractionKernels_hpp__ +#define __grainInteractionKernels_hpp__ + + +#include "grainTriSurfaceContact.hpp" + +namespace pFlow::grainInteractionKernels +{ + +template< + typename ContactForceModel, + typename ContactListType> +struct ppInteractionFunctor +{ + + using PairType = typename ContactListType::PairType; + using ValueType = typename ContactListType::ValueType; + + real dt_; + + ContactForceModel forceModel_; + ContactListType tobeFilled_; + + deviceViewType1D diam_; + deviceViewType1D coarseGrainFactor_; + deviceViewType1D propId_; + deviceViewType1D pos_; + deviceViewType1D lVel_; + deviceViewType1D rVel_; + deviceViewType1D cForce_; + deviceViewType1D cTorque_; + + + ppInteractionFunctor( + real dt, + ContactForceModel forceModel, + ContactListType tobeFilled, + deviceViewType1D diam, + deviceViewType1D coarseGrainFactor, + deviceViewType1D propId, + deviceViewType1D pos, + deviceViewType1D lVel, + deviceViewType1D rVel, + deviceViewType1D cForce, + deviceViewType1D cTorque ) + : + dt_(dt), + forceModel_(forceModel), + tobeFilled_(tobeFilled), + diam_(diam), + coarseGrainFactor_(coarseGrainFactor), + propId_(propId), + pos_(pos), + lVel_(lVel), + rVel_(rVel), + cForce_(cForce), // this is converted to an atomic vector + cTorque_(cTorque) // this is converted to an atomic vector + {} + + INLINE_FUNCTION_HD + void operator()(const uint32 n)const + { + + if(!tobeFilled_.isValid(n))return; + + auto [i,j] = tobeFilled_.getPair(n); + + real Ri = 0.5*diam_[i]; + real Rj = 0.5*diam_[j]; + + real cGFi = coarseGrainFactor_[j]; + real cGFj = coarseGrainFactor_[j]; + + realx3 xi = pos_[i]; + realx3 xj = pos_[j]; + real dist = length(xj-xi); + real ovrlp = (Ri+Rj) - dist; + + if( ovrlp >0.0 ) + { + + auto Vi = lVel_[i]; + auto Vj = lVel_[j]; + auto wi = rVel_[i]; + auto wj = rVel_[j]; + auto Nij = (xj-xi)/dist; + auto Vr = Vi - Vj + cross((Ri*wi+Rj*wj), Nij); + + auto history = tobeFilled_.getValue(n); + + int32 propId_i = propId_[i]; + int32 propId_j = propId_[j]; + + realx3 FCn, FCt, Mri, Mrj, Mij, Mji; + + // calculates contact force + forceModel_.contactForce( + dt_, i, j, + propId_i, propId_j, + Ri, Rj, cGFi , cGFj , + ovrlp, + Vr, Nij, + history, + FCn, FCt + ); + + forceModel_.rollingFriction( + dt_, i, j, + propId_i, propId_j, + Ri, Rj, cGFi , cGFj , + 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(&cForce_[i].x_,FC.x_); + Kokkos::atomic_add(&cForce_[i].y_,FC.y_); + Kokkos::atomic_add(&cForce_[i].z_,FC.z_); + + Kokkos::atomic_add(&cForce_[j].x_,-FC.x_); + Kokkos::atomic_add(&cForce_[j].y_,-FC.y_); + Kokkos::atomic_add(&cForce_[j].z_,-FC.z_); + + Kokkos::atomic_add(&cTorque_[i].x_, Mij.x_); + Kokkos::atomic_add(&cTorque_[i].y_, Mij.y_); + Kokkos::atomic_add(&cTorque_[i].z_, Mij.z_); + + Kokkos::atomic_add(&cTorque_[j].x_, Mji.x_); + Kokkos::atomic_add(&cTorque_[j].y_, Mji.y_); + Kokkos::atomic_add(&cTorque_[j].z_, Mji.z_); + + + tobeFilled_.setValue(n,history); + + } + else + { + tobeFilled_.setValue(n, ValueType()); + } + + } +}; + + +template< + typename ContactForceModel, + typename ContactListType, + typename TraingleAccessor, + typename MotionModel> +struct pwInteractionFunctor +{ + using PairType = typename ContactListType::PairType; + using ValueType = typename ContactListType::ValueType; + + real dt_; + + ContactForceModel forceModel_; + ContactListType tobeFilled_; + + TraingleAccessor triangles_; + MotionModel motionModel_; + + deviceViewType1D diam_; + deviceViewType1D coarseGrainFactor_; + deviceViewType1D propId_; + deviceViewType1D pos_; + deviceViewType1D lVel_; + deviceViewType1D rVel_; + deviceViewType1D cForce_; + deviceViewType1D cTorque_; + deviceViewType1D wTriMotionIndex_; + deviceViewType1D wPropId_; + deviceViewType1D wCForce_; + + + pwInteractionFunctor( + real dt, + ContactForceModel forceModel, + ContactListType tobeFilled, + TraingleAccessor triangles, + MotionModel motionModel , + deviceViewType1D diam , + deviceViewType1D coarseGrainFactor, + deviceViewType1D propId, + deviceViewType1D pos , + deviceViewType1D lVel, + deviceViewType1D rVel, + deviceViewType1D cForce, + deviceViewType1D cTorque , + deviceViewType1D wTriMotionIndex, + deviceViewType1D wPropId, + deviceViewType1D wCForce) + : + dt_(dt), + forceModel_(forceModel), + tobeFilled_(tobeFilled), + triangles_(triangles) , + motionModel_(motionModel) , + diam_(diam) , + coarseGrainFactor_(coarseGrainFactor), + propId_(propId), + pos_(pos) , + lVel_(lVel), + rVel_(rVel) , + cForce_(cForce), + cTorque_(cTorque) , + wTriMotionIndex_(wTriMotionIndex) , + wPropId_(wPropId), + wCForce_(wCForce) + {} + + INLINE_FUNCTION_HD + void operator()(const int32 n)const + { + + if(!tobeFilled_.isValid(n))return; + + auto [i,tj] = tobeFilled_.getPair(n); + + real Ri = 0.5*diam_[i]; + real Rj = 10000.0; + real cGFi = coarseGrainFactor_[i]; + real cGFj = coarseGrainFactor_[i]; + realx3 xi = pos_[i]; + + realx3x3 tri = triangles_(tj); + real ovrlp; + realx3 Nij, cp; + + if( pFlow::grnTriInteraction::isGrainInContactBothSides( + tri, xi, Ri, ovrlp, Nij, cp) ) + { + auto Vi = lVel_[i]; + auto wi = rVel_[i]; + + int32 mInd = wTriMotionIndex_[tj]; + + auto Vw = motionModel_(mInd, cp); + + //output<< "par-wall index "<< i<<" - "<< tj< 1) + { + Ln = 2; return true; + } + else if (nuW >= 0) + { + if (nu >= 0 && w >= 0) + { + Ln = 0; return true; + } + if (nu > 0 && w < 0) + { + Ln = 1; return true; + } + if (nu < 0 && w > 0) + { + Ln = 3; return true; + } + } + else + { + Ln = 1; return true; + } + + return false; +} + +INLINE_FUNCTION_HD +bool isGrainInContactActiveSide( + const realx3& p1, + const realx3& p2, + const realx3& p3, + const realx3& cntr, + real rad, + real& ovrlp, + realx3& norm, + realx3& cp) +{ + + triWall wall(true, p1,p2,p3); + + real dist = wall.normalDistFromWall(cntr); + + if(dist < 0.0 )return false; + + ovrlp = rad - dist; + + if (ovrlp > 0) + { + realx3 ptOnPlane = wall.nearestPointOnWall(cntr); + + if (pointInPlane(p1, p2, p3, ptOnPlane)) + { + cp = ptOnPlane; + norm = -wall.n_; + return true; + } + + realx3 lnv; + + if (pLine(p1,p2).lineGrainCheck(cntr, rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(p2,p3).lineGrainCheck(cntr, rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(p3,p1).lineGrainCheck(cntr, rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + } + + return false; +} + +INLINE_FUNCTION_HD +bool isGrainInContactBothSides( + const realx3x3& tri, + const realx3& cntr, + real Rad, + real& ovrlp, + realx3& norm, + realx3& cp) +{ + + triWall wall(true, tri.x_,tri.y_,tri.z_); + + real dist = wall.normalDistFromWall(cntr); + + + ovrlp = Rad - abs(dist); + + if (ovrlp > 0) + { + realx3 ptOnPlane = wall.nearestPointOnWall(cntr); + + if (pointInPlane(tri.x_,tri.y_,tri.z_,ptOnPlane)) + { + cp = ptOnPlane; + + if(dist >= 0.0) + norm = -wall.n_; + else + norm = wall.n_; + return true; + } + + realx3 lnv; + + if (pLine(tri.x_, tri.y_).lineGrainCheck(cntr, Rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(tri.y_, tri.z_).lineGrainCheck(cntr, Rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + + if ( pLine(tri.z_, tri.x_).lineGrainCheck(cntr, Rad, lnv, cp, ovrlp)) + { + norm = -lnv; + return true; + } + } + + return false; +} + + +} // pFlow::grnTriInteraction + + +#endif //__grainTriSurfaceContact_hpp__ \ No newline at end of file diff --git a/src/Interaction/grainInteraction/grainInteraction/pLine.hpp b/src/Interaction/grainInteraction/grainInteraction/pLine.hpp new file mode 100644 index 00000000..848ff66f --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteraction/pLine.hpp @@ -0,0 +1,107 @@ +/*------------------------------- 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 __pLine_hpp__ +#define __pLine_hpp__ + +#include "types.hpp" + +namespace pFlow::grnTriInteraction +{ + +struct pLine +{ + + realx3 p1_; // point 1 + realx3 p2_; // piont 2 + realx3 v_; // direction vector + real L_; // line lenght + + INLINE_FUNCTION_HD + pLine(){} + + INLINE_FUNCTION_HD + pLine(const realx3 &p1, const realx3 &p2) + : + p1_(p1), + p2_(p2), + v_(p2-p1), + L_(length(v_)) + {} + + // get a point on the line based on input 0<= t <= 1 + INLINE_FUNCTION_HD + realx3 point(real t)const { + return v_ * t + p1_; + } + + // return the projected point of point p on line + INLINE_FUNCTION_HD + realx3 projectPoint(const realx3 &p) const + { + return point(projectNormLength(p)); + } + + // calculates the normalized distance between projected p and p1 + INLINE_FUNCTION_HD + real projectNormLength(realx3 p) const + { + realx3 w = p - p1_; + return dot(w,v_) / (L_*L_); + } + + INLINE_FUNCTION_HD + bool lineGrainCheck( + const realx3 pos, + real Rad, + realx3 &nv, + realx3 &cp, + real &ovrlp)const + { + + + real t = projectNormLength(pos); + + if(t >= 0.0 && t <= 1.0) cp = point(t); + else if(t >= (-Rad / L_) && t < 0.0) cp = point(0.0); + else if(t>1.0 && t >= (1.0 + Rad / L_)) cp = point(1.0); + else return false; + + realx3 vec = pos - cp; // from cp to pos + + real dist = length(vec); + ovrlp = Rad - dist; + + if (ovrlp >= 0.0) + { + if (dist > 0) + nv = vec / dist; + else + nv = v_; + return true; + } + + return false; + } +}; + +} //pFlow::grnTriInteractio + +#endif diff --git a/src/Interaction/grainInteraction/grainInteraction/triWall.hpp b/src/Interaction/grainInteraction/grainInteraction/triWall.hpp new file mode 100644 index 00000000..5a0086a4 --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteraction/triWall.hpp @@ -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. + +-----------------------------------------------------------------------------*/ + +#ifndef __triWall_hpp__ +#define __triWall_hpp__ + +#include "types.hpp" + +namespace pFlow::grnTriInteraction +{ + +struct triWall +{ + + realx3 n_; + real offset_; + + INLINE_FUNCTION_H + triWall(const realx3& p1, const realx3& p2, const realx3& p3) + { + if(!makeWall(p1,p2,p3, n_, offset_)) + { + fatalErrorInFunction<< + "bad input for the wall.\n"; + fatalExit; + } + } + + INLINE_FUNCTION_HD + triWall(bool, const realx3& p1, const realx3& p2, const realx3& p3) + { + makeWall(p1,p2,p3,n_,offset_); + } + + INLINE_FUNCTION_HD + triWall(const realx3x3& tri) + { + makeWall(tri.x_, tri.y_, tri.z_, n_, offset_); + } + + INLINE_FUNCTION_HD + triWall(const triWall&) = default; + + INLINE_FUNCTION_HD + triWall& operator=(const triWall&) = default; + + INLINE_FUNCTION_HD + triWall(triWall&&) = default; + + INLINE_FUNCTION_HD + triWall& operator=(triWall&&) = default; + + INLINE_FUNCTION_HD + ~triWall()=default; + + + INLINE_FUNCTION_HD + real normalDistFromWall(const realx3 &p) const + { + return dot(n_, p) + offset_; + } + + INLINE_FUNCTION_HD + realx3 nearestPointOnWall(const realx3 &p) const + { + real t = -(dot(n_, p) + offset_); + return realx3(n_.x_*t + p.x_, n_.y_*t + p.y_, n_.z_*t + p.z_); + } + + INLINE_FUNCTION_HD static + bool makeWall( + const realx3& p1, + const realx3& p2, + const realx3& p3, + realx3& n, real& offset) + { + n = cross(p2 - p1, p3 - p1); + real len = length(n); + + if (len < 0.00000000001) return false; + n /= len; + offset = -dot(n, p1); + return true; + } + + +}; + +} + +#endif diff --git a/src/Interaction/grainInteraction/grainInteractionsLinearModels.cpp b/src/Interaction/grainInteraction/grainInteractionsLinearModels.cpp new file mode 100755 index 00000000..6fa77730 --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteractionsLinearModels.cpp @@ -0,0 +1,62 @@ +/*------------------------------- 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 "grainInteraction.hpp" +#include "geometryMotions.hpp" +#include "grainContactForceModels.hpp" +#include "unsortedContactList.hpp" +#include "sortedContactList.hpp" +#include "createBoundaryGrainInteraction.hpp" + + + +#define createInteraction(ForceModel,GeomModel) \ + \ + template class pFlow::grainInteraction< \ + ForceModel, \ + GeomModel, \ + pFlow::unsortedContactList>; \ + \ + template class pFlow::grainInteraction< \ + ForceModel, \ + GeomModel, \ + pFlow::sortedContactList>; \ + createBoundaryGrainInteraction(ForceModel, GeomModel) + + +createInteraction(pFlow::cfModels::grainRolling>, pFlow::rotationAxisMotionGeometry); +createInteraction(pFlow::cfModels::grainRolling>, pFlow::rotationAxisMotionGeometry); + +createInteraction(pFlow::cfModels::grainRolling>, pFlow::rotationAxisMotionGeometry); +createInteraction(pFlow::cfModels::grainRolling>, pFlow::rotationAxisMotionGeometry); + + +createInteraction(pFlow::cfModels::grainRolling>, pFlow::stationaryGeometry); +createInteraction(pFlow::cfModels::grainRolling>, pFlow::stationaryGeometry); + +createInteraction(pFlow::cfModels::grainRolling>, pFlow::stationaryGeometry); +createInteraction(pFlow::cfModels::grainRolling>, pFlow::stationaryGeometry); + + +createInteraction(pFlow::cfModels::grainRolling>, pFlow::vibratingMotionGeometry); +createInteraction(pFlow::cfModels::grainRolling>, pFlow::vibratingMotionGeometry); + +createInteraction(pFlow::cfModels::grainRolling>, pFlow::vibratingMotionGeometry); +createInteraction(pFlow::cfModels::grainRolling>, pFlow::vibratingMotionGeometry); \ No newline at end of file diff --git a/src/Interaction/grainInteraction/grainInteractionsNonLinearModModels.cpp b/src/Interaction/grainInteraction/grainInteractionsNonLinearModModels.cpp new file mode 100755 index 00000000..462b5eb3 --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteractionsNonLinearModModels.cpp @@ -0,0 +1,43 @@ +/*------------------------------- 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 "grainInteraction.hpp" +#include "geometryMotions.hpp" +#include "grainContactForceModels.hpp" +#include "unsortedContactList.hpp" +#include "sortedContactList.hpp" +#include "createBoundaryGrainInteraction.hpp" + + + +#define createInteraction(ForceModel,GeomModel) \ + \ + template class pFlow::grainInteraction< \ + ForceModel, \ + GeomModel, \ + pFlow::unsortedContactList>; \ + \ + template class pFlow::grainInteraction< \ + ForceModel, \ + GeomModel, \ + pFlow::sortedContactList>; \ + createBoundaryGrainInteraction(ForceModel, GeomModel) + + diff --git a/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp b/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp new file mode 100755 index 00000000..462b5eb3 --- /dev/null +++ b/src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp @@ -0,0 +1,43 @@ +/*------------------------------- 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 "grainInteraction.hpp" +#include "geometryMotions.hpp" +#include "grainContactForceModels.hpp" +#include "unsortedContactList.hpp" +#include "sortedContactList.hpp" +#include "createBoundaryGrainInteraction.hpp" + + + +#define createInteraction(ForceModel,GeomModel) \ + \ + template class pFlow::grainInteraction< \ + ForceModel, \ + GeomModel, \ + pFlow::unsortedContactList>; \ + \ + template class pFlow::grainInteraction< \ + ForceModel, \ + GeomModel, \ + pFlow::sortedContactList>; \ + createBoundaryGrainInteraction(ForceModel, GeomModel) + + diff --git a/src/Particles/CMakeLists.txt b/src/Particles/CMakeLists.txt index b7f5ddc6..7cc91295 100644 --- a/src/Particles/CMakeLists.txt +++ b/src/Particles/CMakeLists.txt @@ -9,8 +9,17 @@ particles/regularParticleIdHandler/regularParticleIdHandler.cpp SphereParticles/sphereShape/sphereShape.cpp SphereParticles/sphereParticles/sphereParticles.cpp SphereParticles/sphereParticles/sphereParticlesKernels.cpp + +GrainParticles/grainShape/grainShape.cpp +GrainParticles/grainParticles/grainParticles.cpp +GrainParticles/grainParticles/grainParticlesKernels.cpp + SphereParticles/boundarySphereParticles.cpp SphereParticles/boundarySphereParticlesList.cpp + +GrainParticles/boundaryGrainParticles.cpp +GrainParticles/boundaryGrainParticlesList.cpp + Insertion/collisionCheck/collisionCheck.cpp Insertion/insertionRegion/insertionRegion.cpp Insertion/insertion/insertion.cpp diff --git a/src/Particles/GrainParticles/boundaryGrainParticles.cpp b/src/Particles/GrainParticles/boundaryGrainParticles.cpp new file mode 100644 index 00000000..49c8b7b2 --- /dev/null +++ b/src/Particles/GrainParticles/boundaryGrainParticles.cpp @@ -0,0 +1,64 @@ +#include "boundaryGrainParticles.hpp" +#include "boundaryBase.hpp" +#include "grainParticles.hpp" + + +pFlow::boundaryGrainParticles::boundaryGrainParticles( + const boundaryBase &boundary, + grainParticles &prtcls +) +: + generalBoundary(boundary, prtcls.pStruct(), "", ""), + particles_(prtcls) +{ + +} + +pFlow::grainParticles &pFlow::boundaryGrainParticles::Particles() +{ + return particles_; +} + +const pFlow::grainParticles &pFlow::boundaryGrainParticles::Particles() const +{ + return particles_; +} + +pFlow::uniquePtr pFlow::boundaryGrainParticles::create( + const boundaryBase &boundary, + grainParticles &prtcls +) +{ + + word bType = angleBracketsNames2( + "boundaryGrainParticles", + pFlowProcessors().localRunTypeName(), + boundary.type()); + + word altBType{"boundaryGrainParticles"}; + + if( boundaryBasevCtorSelector_.search(bType) ) + { + pOutput.space(4)<<"Creating boundary "<< Green_Text(bType)<< + " for "<"); + + boundaryGrainParticles( + const boundaryBase &boundary, + grainParticles& prtcls + ); + + create_vCtor( + boundaryGrainParticles, + boundaryBase, + ( + const boundaryBase &boundary, + grainParticles& prtcls + ), + (boundary, prtcls) + ); + + add_vCtor( + boundaryGrainParticles, + boundaryGrainParticles, + boundaryBase + ); + + grainParticles& Particles(); + + const grainParticles& Particles()const; + + bool hearChanges( + real t, + real dt, + uint32 iter, + const message &msg, + const anyList &varList) override + { + return true; + } + + virtual + bool acceleration(const timeInfo& ti, const realx3& g) + { + return true; + } + + static + uniquePtr create( + const boundaryBase &boundary, + grainParticles& prtcls); + +}; + + +} + + + +#endif diff --git a/src/Particles/GrainParticles/boundaryGrainParticlesList.cpp b/src/Particles/GrainParticles/boundaryGrainParticlesList.cpp new file mode 100644 index 00000000..4cd941bd --- /dev/null +++ b/src/Particles/GrainParticles/boundaryGrainParticlesList.cpp @@ -0,0 +1,19 @@ +#include "boundaryGrainParticlesList.hpp" + +pFlow::boundaryGrainParticlesList::boundaryGrainParticlesList( + const boundaryList &bndrs, + grainParticles &prtcls +) +: + ListPtr(bndrs.size()), + boundaries_(bndrs) +{ + for(auto i=0; iset + ( + i, + boundaryGrainParticles::create(boundaries_[i], prtcls) + ); + } +} \ No newline at end of file diff --git a/src/Particles/GrainParticles/boundaryGrainParticlesList.hpp b/src/Particles/GrainParticles/boundaryGrainParticlesList.hpp new file mode 100644 index 00000000..4fe22ad8 --- /dev/null +++ b/src/Particles/GrainParticles/boundaryGrainParticlesList.hpp @@ -0,0 +1,36 @@ + + +#ifndef __boundaryGrainParticlesList_hpp__ +#define __boundaryGrainParticlesList_hpp__ + +#include "ListPtr.hpp" +#include "boundaryList.hpp" +#include "boundaryGrainParticles.hpp" + +namespace pFlow +{ + +class boundaryGrainParticlesList +: + public ListPtr +{ +private: + + const boundaryList& boundaries_; + +public: + + boundaryGrainParticlesList( + const boundaryList& bndrs, + grainParticles& prtcls + ); + + ~boundaryGrainParticlesList()=default; + +}; + +} + + + +#endif \ No newline at end of file diff --git a/src/Particles/GrainParticles/grainParticles/grainParticles.cpp b/src/Particles/GrainParticles/grainParticles/grainParticles.cpp new file mode 100644 index 00000000..8293a510 --- /dev/null +++ b/src/Particles/GrainParticles/grainParticles/grainParticles.cpp @@ -0,0 +1,422 @@ +/*------------------------------- 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 "grainParticles.hpp" +#include "systemControl.hpp" +#include "vocabs.hpp" +#include "grainParticlesKernels.hpp" + +bool pFlow::grainParticles::initializeParticles() +{ + + using exeSpace = typename realPointField_D::execution_space; + using policy = Kokkos::RangePolicy< + exeSpace, + Kokkos::IndexType>; + + auto [minIndex, maxIndex] = minMax(shapeIndex().internal()); + + if( !grains_.indexValid(maxIndex) ) + { + fatalErrorInFunction<< + "the maximum value of shapeIndex is "<< maxIndex << + " which is not valid."<(shapeNames.size()); + + propIds.clear(); + propIds.reserve(numNew); + + diams.clear(); + diams.reserve(numNew); + + coarseGrainFactors.clear(); + coarseGrainFactors.reserve(numNew); + + m.clear(); + m.reserve(numNew); + + Is.clear(); + Is.reserve(numNew); + + shIndex.clear(); + shIndex.reserve(numNew); + + + for(const auto& name:shapeNames) + { + uint32 indx; + if(grains_.shapeNameToIndex(name,indx)) + { + shIndex.push_back(indx); + Is.push_back( grains_.Inertia(indx)); + m.push_back(grains_.mass(indx)); + diams.push_back(grains_.boundingDiameter(indx)); + coarseGrainFactors.push_back(grains_.coarseGrainFactor(indx)); + propIds.push_back( grains_.propertyId(indx)); + } + else + { + fatalErrorInFunction<<"Shape name "<< name << + "does not exist. The list is "<(0.0000000001) + ), + rVelocity_ + ( + objectFile + ( + "rVelocity", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + dynPointStruct(), + zero3 + ), + rAcceleration_ + ( + objectFile( + "rAcceleration", + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS + ), + dynPointStruct(), + zero3 + ), + boundaryGrainParticles_ + ( + dynPointStruct().boundaries(), + *this + ), + accelerationTimer_( + "Acceleration", &this->timers() ), + intPredictTimer_( + "Integration-predict", &this->timers() ), + intCorrectTimer_( + "Integration-correct", &this->timers() ), + fieldUpdateTimer_( + "fieldUpdate", &this->timers() ) +{ + + auto intMethod = control.settingsDict().getVal("integrationMethod"); + REPORT(1)<<"Creating integration method "<dt(); + dynPointStruct().predict(dt, accelertion()); + rVelIntegration_().predict(dt,rVelocity_, rAcceleration_); + intPredictTimer_.end(); + + fieldUpdateTimer_.start(); + propertyId_.updateBoundariesSlaveToMasterIfRequested(); + grainDiameter_.updateBoundariesSlaveToMasterIfRequested(); + coarseGrainFactor_.updateBoundariesSlaveToMasterIfRequested(); + mass_.updateBoundariesSlaveToMasterIfRequested(); + I_.updateBoundariesSlaveToMasterIfRequested(); + rVelocity_.updateBoundariesSlaveToMasterIfRequested(); + rAcceleration_.updateBoundariesSlaveToMasterIfRequested(); + rVelIntegration_().updateBoundariesSlaveToMasterIfRequested(); + fieldUpdateTimer_.end(); + + return true; +} + +bool pFlow::grainParticles::iterate() +{ + + timeInfo ti = TimeInfo(); + realx3 g = control().g(); + + particles::iterate(); + accelerationTimer_.start(); + pFlow::grainParticlesKernels::acceleration( + g, + mass().deviceViewAll(), + contactForce().deviceViewAll(), + I().deviceViewAll(), + contactTorque().deviceViewAll(), + dynPointStruct().activePointsMaskDevice(), + accelertion().deviceViewAll(), + rAcceleration().deviceViewAll() + ); + for(auto& bndry:boundaryGrainParticles_) + { + bndry->acceleration(ti, g); + } + accelerationTimer_.end(); + + intCorrectTimer_.start(); + + if(!dynPointStruct().correct(dt(), accelertion())) + { + return false; + } + if(!rVelIntegration_().correct( + dt(), + rVelocity_, + rAcceleration_)) + { + return false; + } + + intCorrectTimer_.end(); + + return true; +} + +bool pFlow::grainParticles::insertParticles +( + const realx3Vector &position, + const wordVector &shapesNames, + const anyList &setVarList +) +{ + anyList newVarList(setVarList); + + realVector mass("mass"); + realVector I("I"); + realVector diameter("diameter"); + realVector coarseGrainFactor("coarseGrainFactor"); + uint32Vector propId("propId"); + uint32Vector shapeIdx("shapeIdx"); + + if(!getParticlesInfoFromShape( + shapesNames, + propId, + diameter, + coarseGrainFactor, + mass, + I, + shapeIdx)) + { + return false; + } + + newVarList.emplaceBack( + mass_.name()+"Vector", + std::move(mass)); + + newVarList.emplaceBack( + I_.name()+"Vector", + std::move(I)); + + newVarList.emplaceBack( + grainDiameter_.name()+"Vector", + std::move(diameter)); + + newVarList.emplaceBack( + coarseGrainFactor_.name()+"Vector", + std::move(coarseGrainFactor)); + + newVarList.emplaceBack( + propertyId_.name()+"Vector", + std::move(propId)); + + newVarList.emplaceBack( + shapeIndex().name()+"Vector", + std::move(shapeIdx)); + + if(!dynPointStruct().insertPoints(position, newVarList)) + { + return false; + } + + return true; +} + +pFlow::word pFlow::grainParticles::shapeTypeName()const +{ + return "grain"; +} + +const pFlow::shape &pFlow::grainParticles::getShapes() const +{ + return grains_; +} + +void pFlow::grainParticles::boundingSphereMinMax +( + real & minDiam, + real& maxDiam +)const +{ + minDiam = grains_.minBoundingSphere(); + maxDiam = grains_.maxBoundingSphere(); +} diff --git a/src/Particles/GrainParticles/grainParticles/grainParticles.hpp b/src/Particles/GrainParticles/grainParticles/grainParticles.hpp new file mode 100644 index 00000000..9a0a35ac --- /dev/null +++ b/src/Particles/GrainParticles/grainParticles/grainParticles.hpp @@ -0,0 +1,236 @@ +/*------------------------------- 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. + +-----------------------------------------------------------------------------*/ +/** + * @class pFlow::sphereParticles + * + * @brief Class for managing spherical particles + * + * This is a top-level class that contains the essential components for + * defining spherical prticles in a DEM simulation. + */ + +#ifndef __grainParticles_hpp__ +#define __grainParticles_hpp__ + +#include "indexContainer.hpp" +#include "particles.hpp" +#include "property.hpp" +#include "grainShape.hpp" +#include "boundaryGrainParticlesList.hpp" +#include "systemControl.hpp" + +namespace pFlow +{ + +class grainParticles : public particles +{ +public: + + using ShapeType = grainShape; + +private: + + /// reference to shapes + ShapeType grains_; + + /// property id on device + uint32PointField_D propertyId_; + + /// diameter / boundig sphere size of particles on device + realPointField_D grainDiameter_; + + realPointField_D coarseGrainFactor_; + + + /// mass of particles field + realPointField_D mass_; + + /// pointField of inertial of particles + realPointField_D I_; + + /// pointField of rotational Velocity of particles on device + realx3PointField_D rVelocity_; + + /// pointField of rotational acceleration of particles on device + realx3PointField_D rAcceleration_; + + /// boundaries + boundaryGrainParticlesList boundaryGrainParticles_; + + /// rotational velocity integrator + uniquePtr rVelIntegration_ = nullptr; + + /// timer for acceleration computations + Timer accelerationTimer_; + + /// timer for integration computations (prediction step) + Timer intPredictTimer_; + + /// timer for integration computations (correction step) + Timer intCorrectTimer_; + + Timer fieldUpdateTimer_; + +private: + + + bool getParticlesInfoFromShape( + const wordVector& shapeNames, + uint32Vector& propIds, + realVector& diams, + realVector& coarseGrainFactors, + realVector& m, + realVector& Is, + uint32Vector& shIndex + ); + /*bool initializeParticles(); + + bool insertSphereParticles( + const wordVector& names, + const int32IndexContainer& indices, + bool setId = true); + + virtual uniquePtr> getFieldObjectList()const override; + */ + +public: + + /// construct from systemControl and property + grainParticles(systemControl& control, const property& prop); + + ~grainParticles() override = default; + + /** + * Insert new particles in position with specified shapes + * + * This function is involked by inserted object to insert new set of + * particles into the simulation. \param position position of new particles + * \param shape shape of new particles + * \param setField initial value of the selected fields for new particles + */ + /*bool insertParticles + ( + const realx3Vector& position, + const wordVector& shapes, + const setFieldList& setField + ) override ;*/ + + // TODO: make this method private later + bool initializeParticles(); + + /// const reference to shapes object + const auto& grains() const + { + return grains_; + } + + /// const reference to inertia pointField + const auto& I() const + { + return I_; + } + + /// reference to inertia pointField + auto& I() + { + return I_; + } + + const auto& rVelocity() const + { + return rVelocity_; + } + + auto& rVelocity() + { + return rVelocity_; + } + + bool hearChanges( + real t, + real dt, + uint32 iter, + const message& msg, + const anyList& varList + ) override + { + notImplementedFunction; + return false; + } + + const uint32PointField_D& propertyId() const override + { + return propertyId_; + } + + const realPointField_D& diameter() const override + { + return grainDiameter_; + } + + const realPointField_D& coarseGrainFactor() const + { + return coarseGrainFactor_; + } + + + const realPointField_D& mass() const override + { + return mass_; + } + + /// before iteration step + bool beforeIteration() override; + + /// iterate particles + bool iterate() override; + + bool insertParticles( + const realx3Vector& position, + const wordVector& shapesNames, + const anyList& setVarList + ) override; + + realx3PointField_D& rAcceleration() override + { + return rAcceleration_; + } + + const realx3PointField_D& rAcceleration() const override + { + return rAcceleration_; + } + + const realPointField_D& boundingSphere() const override + { + return diameter(); + } + + word shapeTypeName() const override; + + const shape& getShapes() const override; + + void boundingSphereMinMax(real& minDiam, real& maxDiam) const override; + +};// grainParticles + +} // pFlow + +#endif //__sphereParticles_hpp__ diff --git a/src/Particles/GrainParticles/grainParticles/grainParticlesKernels.cpp b/src/Particles/GrainParticles/grainParticles/grainParticlesKernels.cpp new file mode 100644 index 00000000..51891205 --- /dev/null +++ b/src/Particles/GrainParticles/grainParticles/grainParticlesKernels.cpp @@ -0,0 +1,101 @@ +/*------------------------------- 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 "grainParticlesKernels.hpp" + +using policy = Kokkos::RangePolicy< + pFlow::DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType>; + +void pFlow::grainParticlesKernels::addMassDiamInertiaProp +( + deviceViewType1D shapeIndex, + deviceViewType1D mass, + deviceViewType1D diameter, + deviceViewType1D coarseGrainFactor, + deviceViewType1D I, + deviceViewType1D propertyId, + pFlagTypeDevice incld, + deviceViewType1D src_mass, + deviceViewType1D src_diameter, + deviceViewType1D src_I, + deviceViewType1D src_propertyId +) +{ + auto aRange = incld.activeRange(); + + Kokkos::parallel_for( + "particles::initInertia", + policy(aRange.start(), aRange.end()), + LAMBDA_HD(uint32 i) + { + if(incld(i)) + { + uint32 index = shapeIndex[i]; + I[i] = src_I[index]; + diameter[i] = src_diameter[index]; + mass[i] = src_mass[index]; + propertyId[i] = src_propertyId[index]; + } + }); + +} + +void pFlow::grainParticlesKernels::acceleration +( + const realx3& g, + const deviceViewType1D& mass, + const deviceViewType1D& force, + const deviceViewType1D& I, + const deviceViewType1D& torque, + const pFlagTypeDevice& incld, + deviceViewType1D lAcc, + deviceViewType1D rAcc +) +{ + + auto activeRange = incld.activeRange(); + if(incld.isAllActive()) + { + Kokkos::parallel_for( + "pFlow::grainParticlesKernels::acceleration", + policy(activeRange.start(), activeRange.end()), + LAMBDA_HD(uint32 i){ + lAcc[i] = force[i]/mass[i] + g; + rAcc[i] = torque[i]/I[i]; + }); + } + else + { + Kokkos::parallel_for( + "pFlow::grainParticlesKernels::acceleration", + policy(activeRange.start(), activeRange.end()), + LAMBDA_HD(uint32 i){ + if(incld(i)) + { + lAcc[i] = force[i]/mass[i] + g; + rAcc[i] = torque[i]/I[i]; + } + }); + + } + + Kokkos::fence(); +} diff --git a/src/Particles/GrainParticles/grainParticles/grainParticlesKernels.hpp b/src/Particles/GrainParticles/grainParticles/grainParticlesKernels.hpp new file mode 100644 index 00000000..daed7808 --- /dev/null +++ b/src/Particles/GrainParticles/grainParticles/grainParticlesKernels.hpp @@ -0,0 +1,59 @@ +/*------------------------------- 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 __grainParticlesKernels_hpp__ +#define __grainParticlesKernels_hpp__ + +#include "types.hpp" +#include "pointFlag.hpp" + +namespace pFlow::grainParticlesKernels +{ + +void addMassDiamInertiaProp( + deviceViewType1D shapeIndex, + deviceViewType1D mass, + deviceViewType1D diameter, + deviceViewType1D coarseGrainFactor, + + deviceViewType1D I, + deviceViewType1D propertyId, + pFlagTypeDevice incld, + deviceViewType1D src_mass, + deviceViewType1D src_grainDiameter, + deviceViewType1D src_I, + deviceViewType1D src_propertyId +); + +void acceleration( + const realx3& g, + const deviceViewType1D& mass, + const deviceViewType1D& force, + const deviceViewType1D& I, + const deviceViewType1D& torque, + const pFlagTypeDevice& incld, + deviceViewType1D lAcc, + deviceViewType1D rAcc +); + + +} + +#endif diff --git a/src/Particles/GrainParticles/grainShape/grainShape.cpp b/src/Particles/GrainParticles/grainShape/grainShape.cpp new file mode 100644 index 00000000..4558940b --- /dev/null +++ b/src/Particles/GrainParticles/grainShape/grainShape.cpp @@ -0,0 +1,242 @@ +/*------------------------------- 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 "grainShape.hpp" + + +bool pFlow::grainShape::readFromDictionary3() +{ + + grainDiameters_ = getVal("grainDiameters"); + + sphereDiameters_ = getVal("sphereDiameters"); + + coarseGrainFactor_ = grainDiameters_ / sphereDiameters_ ; + + + if(grainDiameters_.size() != numShapes() ) + { + fatalErrorInFunction<< + " number of elements in grain diameters in "<< globalName()<<" is not consistent"<"); + + grainShape( + const word& fileName, + repository* owner, + const property& prop); + + + ~grainShape() override = default; + + //// - Methods + + real maxBoundingSphere()const override; + + real minBoundingSphere()const override; + + bool boundingDiameter(uint32 index, real& bDiam)const override; + + real boundingDiameter(uint32 index)const override; + + realVector boundingDiameter()const override; + + real coarseGrainFactor(uint32 index)const ; + + realVector coarseGrainFactor()const ; + + real orginalDiameter(uint32 index)const ; + + realVector orginalDiameter()const ; + + bool mass(uint32 index, real& m)const override; + + real mass(uint32 index) const override; + + realVector mass()const override; + + realVector density() const override; + + bool Inertia(uint32 index, real& I)const override; + + real Inertia(uint32 index)const override; + + realVector Inertia()const override; + + bool Inertia_xx(uint32 index, real& Ixx)const override; + + real Inertial_xx(uint32 index)const override; + + bool Inertia_yy(uint32 index, real& Iyy)const override; + + real Inertial_yy(uint32 index)const override; + + bool Inertia_zz(uint32 index, real& Izz)const override; + + real Inertial_zz(uint32 index)const override; + +}; + +} // pFlow + +#endif //__grainShape_hpp__ diff --git a/src/Particles/Insertion/Insertion/Insertions.cpp b/src/Particles/Insertion/Insertion/Insertions.cpp index 2bc3bd92..5405b3f6 100644 --- a/src/Particles/Insertion/Insertion/Insertions.cpp +++ b/src/Particles/Insertion/Insertion/Insertions.cpp @@ -21,4 +21,5 @@ Licence: #include "Insertions.hpp" -template class pFlow::Insertion; \ No newline at end of file +template class pFlow::Insertion; +template class pFlow::Insertion; \ No newline at end of file diff --git a/src/Particles/Insertion/Insertion/Insertions.hpp b/src/Particles/Insertion/Insertion/Insertions.hpp index deaafba2..b513295d 100644 --- a/src/Particles/Insertion/Insertion/Insertions.hpp +++ b/src/Particles/Insertion/Insertion/Insertions.hpp @@ -24,11 +24,15 @@ Licence: #include "Insertion.hpp" #include "sphereShape.hpp" +#include "grainShape.hpp" + namespace pFlow { using sphereInsertion = Insertion ; +using grainInsertion = Insertion ; + } diff --git a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp index 9a3d13c0..86e71b89 100644 --- a/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp +++ b/src/phasicFlow/structuredData/boundaries/boundaryBase/boundaryBase.hpp @@ -164,10 +164,12 @@ protected: return true; } - uint32 markInNegativeSide(const word& name, uint32Vector_D& markedIndices )const; public: + uint32 markInNegativeSide(const word& name, uint32Vector_D& markedIndices )const; + + TypeInfo("boundaryBase"); boundaryBase(