xfygogo/phasicFlow
1
0
Fork 0
mirror of https://github.com/PhasicFlow/phasicFlow.git synced 2025-06-22 16:28:30 +00:00
Code Issues Packages Projects Releases Wiki Activity

AB5 is added and bug is resolved when particles are being inserted

Browse Source
This commit is contained in:
HRN
2025-02-08 18:06:30 +03:30
parent b9ab015eb1
commit db9b1e62e4
13 changed files with 244 additions and 247 deletions
Download Patch File Download Diff File Expand all files Collapse all files

239
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
(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
)
{
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(
bool intAllActive(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
range activeRng)
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.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;
}