xfygogo
/
phasicFlow
mirror of https://github.com/PhasicFlow/phasicFlow.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #161 from PhasicFlow/develop 

AB5 is added and bug is resolved when particles are being inserted
This commit is contained in:
PhasicFlow 2025-02-08 18:07:39 +03:30 committed by GitHub
parent 3aff0f1fc6 db9b1e62e4
commit cb6d567dab
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
13 changed files with 244 additions and 247 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/Integration/AdamsBashforth2/AdamsBashforth2.cpp
View File

@ -50,7 +50,7 @@ bool intAllActive(
"AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] = damping*(d_y[i] + dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]));
d_y[i] += damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i];
});
Kokkos::fence();
@ -80,7 +80,7 @@ bool intScattered
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] = damping*(d_y[i] + dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]));
d_y[i] += damping * dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_dy1[i] = d_dy[i];
}
});
@ -155,11 +155,11 @@ bool pFlow::AdamsBashforth2::correct
bool success = false;
if(dy1l.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1l, damping);
success = intAllActive(dt, y.field(), dy, dy1(), damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1l, damping);
success = intScattered(dt, y.field(), dy, dy1(), damping);
}
success = success && boundaryList_.correct(dt, y, dy);
@ -177,11 +177,11 @@ bool pFlow::AdamsBashforth2::correctPStruct(
bool success = false;
if(dy1l.isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1l);
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1l);
success = intScattered(dt, pStruct.pointPosition(), vel, dy1());
}
success = success && boundaryList_.correctPStruct(dt, pStruct, vel);
@ -189,7 +189,7 @@ bool pFlow::AdamsBashforth2::correctPStruct(
return success;
}
bool pFlow::AdamsBashforth2::hearChanges
/*bool pFlow::AdamsBashforth2::hearChanges
(
const timeInfo &ti,
const message &msg,
@ -206,4 +206,4 @@ bool pFlow::AdamsBashforth2::hearChanges
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}
}*/

4
src/Integration/AdamsBashforth2/AdamsBashforth2.hpp
View File

@ -124,12 +124,12 @@ public:
pointStructure& pStruct,
realx3PointField_D& vel) override;
bool hearChanges
/*bool hearChanges
(
const timeInfo& ti,
const message& msg,
const anyList& varList
) override;
) override;*/
};

14
src/Integration/AdamsBashforth3/AdamsBashforth3.cpp
View File

@ -49,9 +49,9 @@ bool intAllActive(
"AdamsBashforth3::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] = damping*( d_y[i]+ dt*( static_cast<real>(23.0 / 12.0) * d_dy[i]
d_y[i] += damping* dt*( static_cast<real>(23.0 / 12.0) * d_dy[i]
- static_cast<real>(16.0 / 12.0) * d_dy1[i]
+ static_cast<real>(5.0 / 12.0) * d_dy2[i]) );
+ static_cast<real>(5.0 / 12.0) * d_dy2[i]) ;
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
@ -81,14 +81,14 @@ bool intScattered
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
"AdamsBashforth3::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] = damping * (d_y[i] + dt*( static_cast<real>(23.0 / 12.0) * d_dy[i]
d_y[i] += damping * dt*( static_cast<real>(23.0 / 12.0) * d_dy[i]
- static_cast<real>(16.0 / 12.0) * d_dy1[i]
+ static_cast<real>(5.0 / 12.0) * d_dy2[i]));
+ static_cast<real>(5.0 / 12.0) * d_dy2[i]);
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
@ -184,7 +184,7 @@ bool pFlow::AdamsBashforth3::correctPStruct(
return success;
}
bool pFlow::AdamsBashforth3::hearChanges
/*bool pFlow::AdamsBashforth3::hearChanges
(
const timeInfo &ti,
const message &msg,
@ -200,4 +200,4 @@ bool pFlow::AdamsBashforth3::hearChanges
{
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}
}*/

4
src/Integration/AdamsBashforth3/AdamsBashforth3.hpp
View File

@ -108,12 +108,12 @@ public:
realx3PointField_D& vel) override;
bool hearChanges
/*bool hearChanges
(
const timeInfo& ti,
const message& msg,
const anyList& varList
) override;
) override;*/
};

25
src/Integration/AdamsBashforth4/AdamsBashforth4.cpp
View File

@ -48,15 +48,15 @@ bool intAllActive(
auto activeRng = dy1.activeRange();
Kokkos::parallel_for(
"AdamsBashforth3::correct",
"AdamsBashforth4::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] = damping *(d_y[i] + dt*(
d_y[i] += damping * dt*(
static_cast<real>(55.0 / 24.0) * d_dy[i]
- static_cast<real>(59.0 / 24.0) * d_dy1[i]
+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
- static_cast<real>( 9.0 / 24.0) * d_dy3[i]
));
- static_cast<real>( 9.0 / 24.0) * d_dy3[i]);
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
@ -87,17 +87,16 @@ bool intScattered
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
"AdamsBashforth4::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] = damping* ( d_y[i] + dt*(
d_y[i] += damping* dt*(
static_cast<real>(55.0 / 24.0) * d_dy[i]
- static_cast<real>(59.0 / 24.0) * d_dy1[i]
+ static_cast<real>(37.0 / 24.0) * d_dy2[i]
- static_cast<real>( 9.0 / 24.0) * d_dy3[i]
));
- static_cast<real>( 9.0 / 24.0) * d_dy3[i] );
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
@ -133,11 +132,7 @@ pFlow::AdamsBashforth4::AdamsBashforth4
zero3,
zero3
)
{
}
{}
void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
{
@ -190,7 +185,7 @@ bool pFlow::AdamsBashforth4::correctPStruct(
return success;
}
bool pFlow::AdamsBashforth4::hearChanges
/*bool pFlow::AdamsBashforth4::hearChanges
(
const timeInfo &ti,
const message &msg,
@ -207,4 +202,4 @@ bool pFlow::AdamsBashforth4::hearChanges
{
return realx3PointField_D::hearChanges(ti, msg, varList);
}
}
}*/

37
src/Integration/AdamsBashforth4/AdamsBashforth4.hpp
View File

@ -104,50 +104,17 @@ public:
pointStructure& pStruct,
realx3PointField_D& vel) override;
bool hearChanges
/*bool hearChanges
(
const timeInfo& ti,
const message& msg,
const anyList& varList
) override;
) override;*/
};
/*template<typename activeFunctor>
bool pFlow::AdamsBashforth4::intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP )
{
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto activeRng = activeP.activeRange();
Kokkos::parallel_for(
"AdamsBashforth4::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
if( activeP(i))
{
d_y[i] += dt*(
static_cast<real>(55.0 / 24.0) * d_dy[i]
- static_cast<real>(59.0 / 24.0) * d_history[i].dy1_
+ static_cast<real>(37.0 / 24.0) * d_history[i].dy2_
- static_cast<real>( 9.0 / 24.0) * d_history[i].dy3_
);
d_history[i].dy3_ = d_history[i].dy2_;
d_history[i].dy2_ = d_history[i].dy1_;
d_history[i].dy1_ = d_dy[i];
}
});
Kokkos::fence();
return true;
}*/
} // pFlow

235
src/Integration/AdamsBashforth5/AdamsBashforth5.cpp
View File

@ -20,93 +20,178 @@ Licence:
#include "AdamsBashforth5.hpp"
pFlow::AdamsBashforth5::AdamsBashforth5
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
)
:
integration(baseName, owner, pStruct, method),
history_(
owner.emplaceObject<pointField<VectorSingle,AB5History>>(
objectFile(
groupNames(baseName,"AB5History"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
AB5History({zero3,zero3, zero3})))
namespace pFlow
{
}
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool pFlow::AdamsBashforth5::predict
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
{
return true;
}
bool pFlow::AdamsBashforth5::correct
(
bool intAllActive(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
)
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
realx3PointField_D& dy4,
real damping = 1.0)
{
if(this->pStruct().allActive())
{
return intAll(dt, y, dy, this->pStruct().activeRange());
}
else
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
}
return true;
}
bool pFlow::AdamsBashforth5::setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y)
{
return true;
}
bool pFlow::AdamsBashforth5::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
{
auto d_dy = dy.deviceViewAll();
auto d_y = y.deviceViewAll();
auto d_history = history_.deviceViewAll();
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1= dy1.deviceView();
auto d_dy2= dy2.deviceView();
auto d_dy3= dy3.deviceView();
auto d_dy4= dy4.deviceView();
auto activeRng = dy1.activeRange();
Kokkos::parallel_for(
"AdamsBashforth5::correct",
rpIntegration (activeRng.first, activeRng.second),
LAMBDA_HD(int32 i){
d_y[i] += dt*(
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] += damping* dt*(
static_cast<real>(1901.0 / 720.0) * d_dy[i]
- static_cast<real>(2774.0 / 720.0) * d_history[i].dy1_
+ static_cast<real>(2616.0 / 720.0) * d_history[i].dy2_
- static_cast<real>(1274.0 / 720.0) * d_history[i].dy3_
+ static_cast<real>( 251.0 / 720.0) * d_history[i].dy4_
);
d_history[i] = {d_dy[i] ,d_history[i].dy1_, d_history[i].dy2_, d_history[i].dy3_};
- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
d_dy4[i] = d_dy3[i];
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
});
Kokkos::fence();
return true;
}
bool intScattered
(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
realx3PointField_D& dy4,
real damping = 1.0
)
{
auto d_dy = dy.deviceView();
auto d_y = y.deviceView();
auto d_dy1 = dy1.deviceView();
auto d_dy2 = dy2.deviceView();
auto d_dy3 = dy3.deviceView();
auto d_dy4 = dy4.deviceView();
auto activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth5::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] += damping* dt*(
static_cast<real>(1901.0 / 720.0) * d_dy[i]
- static_cast<real>(2774.0 / 720.0) * d_dy1[i]
+ static_cast<real>(2616.0 / 720.0) * d_dy2[i]
- static_cast<real>(1274.0 / 720.0) * d_dy3[i]
+ static_cast<real>( 251.0 / 720.0) * d_dy4[i]);
d_dy4[i] = d_dy3[i];
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
pFlow::AdamsBashforth5::AdamsBashforth5
(
const word &baseName,
pointStructure &pStruct,
const word &method,
const realx3Field_D &initialValField
)
:
AdamsBashforth4(baseName, pStruct, method, initialValField),
dy4_
(
objectFile
(
groupNames(baseName,"dy4"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{}
void pFlow::AdamsBashforth5::updateBoundariesSlaveToMasterIfRequested()
{
AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested();
dy4_.updateBoundariesSlaveToMasterIfRequested();
}
bool pFlow::AdamsBashforth5::correct
(
real dt,
realx3PointField_D &y,
realx3PointField_D &dy,
real damping
)
{
bool success = false;
if(y.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), dy4(), damping);
}
success = success && boundaryList().correct(dt, y, dy);
return success;
}
bool pFlow::AdamsBashforth5::correctPStruct
(
real dt,
pointStructure &pStruct,
realx3PointField_D &vel
)
{
bool success = false;
if(dy2().isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3(), dy4());
}
success = success && boundaryList().correctPStruct(dt, pStruct, vel);
return success;
}

130
src/Integration/AdamsBashforth5/AdamsBashforth5.hpp
View File

@ -22,54 +22,12 @@ Licence:
#define __AdamsBashforth5_hpp__
#include "integration.hpp"
#include "AdamsBashforth4.hpp"
#include "pointFields.hpp"
namespace pFlow
{
struct AB5History
{
TypeInfoNV("AB5History");
realx3 dy1_={0,0,0};
realx3 dy2_={0,0,0};
realx3 dy3_={0,0,0};
realx3 dy4_={0,0,0};
};
INLINE_FUNCTION
iIstream& operator>>(iIstream& str, AB5History& ab5)
{
str.readBegin("AB5History");
str >> ab5.dy1_;
str >> ab5.dy2_;
str >> ab5.dy3_;
str >> ab5.dy4_;
str.readEnd("AB5History");
str.check(FUNCTION_NAME);
return str;
}
INLINE_FUNCTION
iOstream& operator<<(iOstream& str, const AB5History& ab5)
{
str << token::BEGIN_LIST << ab5.dy1_
<< token::SPACE << ab5.dy2_
<< token::SPACE << ab5.dy3_
<< token::SPACE << ab5.dy4_
<< token::END_LIST;
str.check(FUNCTION_NAME);
return str;
}
/**
* Fifth order Adams-Bashforth integration method for solving ODE
@ -78,19 +36,25 @@ iOstream& operator<<(iOstream& str, const AB5History& ab5)
*/
class AdamsBashforth5
:
public integration
public AdamsBashforth4
{
private:
friend class processorAB4BoundaryIntegration;
realx3PointField_D dy4_;
protected:
/// Integration history
pointField<VectorSingle,AB5History>& history_;
const auto& dy4()const
{
return dy4_;
}
/// Range policy for integration kernel
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
auto& dy4()
{
return dy4_;
}
public:
@ -99,20 +63,19 @@ public:
// - Constructors
/// Construct from components
AdamsBashforth5(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method);
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField);
uniquePtr<integration> clone()const override
{
return makeUnique<AdamsBashforth5>(*this);
}
virtual ~AdamsBashforth5()=default;
/// Add this to the virtual constructor table
/// Destructor
~AdamsBashforth5() override =default;
/// Add a this to the virtual constructor table
add_vCtor(
integration,
AdamsBashforth5,
@ -121,44 +84,29 @@ public:
// - Methods
bool predict(
real UNUSED(dt),
realx3Vector_D & UNUSED(y),
realx3Vector_D& UNUSED(dy)) override;
void updateBoundariesSlaveToMasterIfRequested()override;
/// return integration method
word method()const override
{
return "AdamsBashforth5";
}
bool correct(
real dt,
realx3Vector_D & y,
realx3Vector_D& dy) override;
realx3PointField_D& y,
realx3PointField_D& dy,
real damping = 1.0) override;
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) override;
bool needSetInitialVals()const override
{
return false;
}
/// Integrate on all points in the active range
bool intAll(
bool correctPStruct(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng);
/// Integrate on active points in the active range
template<typename activeFunctor>
bool intRange(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
activeFunctor activeP );
pointStructure& pStruct,
realx3PointField_D& vel) override;
};
template<typename activeFunctor>
/*template<typename activeFunctor>
bool pFlow::AdamsBashforth5::intRange(
real dt,
realx3Vector_D& y,
@ -190,7 +138,7 @@ bool pFlow::AdamsBashforth5::intRange(
Kokkos::fence();
return true;
}
}*/
} // pFlow

2
src/Integration/CMakeLists.txt
View File

@ -6,7 +6,7 @@ boundaries/boundaryIntegrationList.cpp
AdamsBashforth2/AdamsBashforth2.cpp
AdamsBashforth3/AdamsBashforth3.cpp
AdamsBashforth4/AdamsBashforth4.cpp
#AdamsBashforth5/AdamsBashforth5.cpp
AdamsBashforth5/AdamsBashforth5.cpp
#AdamsMoulton3/AdamsMoulton3.cpp
#AdamsMoulton4/AdamsMoulton4.cpp

6
src/Interaction/contactSearch/ContactSearch/ContactSearch.hpp
View File

@ -60,9 +60,9 @@ private:
bool force = false) override
{
const auto& position = Particles().pointPosition().deviceViewAll();
const auto& flags = Particles().dynPointStruct().activePointsMaskDevice();
const auto& diam = Particles().boundingSphere().deviceViewAll();
auto position = Particles().pointPosition().deviceViewAll();
auto flags = Particles().dynPointStruct().activePointsMaskDevice();
auto diam = Particles().boundingSphere().deviceViewAll();
return ppwContactSearch_().broadSearch(
ti.iter(),

4
src/Interaction/contactSearch/methods/cellBased/NBS/mapperNBSKernels.cpp
View File

@ -22,7 +22,6 @@ Licence:
#include "mapperNBSKernels.hpp"
void pFlow::mapperNBSKernels::findPointExtends
(
const deviceViewType1D<realx3>& points,
@ -154,6 +153,7 @@ bool pFlow::mapperNBSKernels::buildLists
)
{
auto aRange = flags.activeRange();
auto pp = points;
if(flags.isAllActive() )
{
Kokkos::parallel_for
@ -162,7 +162,7 @@ bool pFlow::mapperNBSKernels::buildLists
deviceRPolicyStatic(aRange.start(), aRange.end()),
LAMBDA_HD(uint32 i)
{
auto ind = searchCell.pointIndex(points[i]);
auto ind = searchCell.pointIndex(pp[i]);
uint32 old = Kokkos::atomic_exchange(&head(ind.x(), ind.y(), ind.z()), i);
next[i] = old;
}

1
src/phasicFlow/containers/pointField/internalField/internalField.cpp
View File

@ -21,6 +21,7 @@ Licence:
template<class T, class MemorySpace>
bool pFlow::internalField<T, MemorySpace>::insert(const anyList& varList)
{
const word eventName = message::eventName(message::ITEMS_INSERT);
if(!varList.checkObjectType<uint32IndexContainer>(eventName))

3
src/phasicFlow/structuredData/pointStructure/internalPoints/internalPoints.cpp
View File

@ -304,7 +304,8 @@ pFlow::internalPoints::insertPoints(
{
return false;
}
//output<<" points \n"<<points<<endl;
//output<<pointPosition_<<endl;
if(!notify(
time().TimeInfo(),
msg,