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

Merge pull request #151 from PhasicFlow/main 

update from main branch
This commit is contained in:
PhasicFlow 2025-01-24 21:12:53 +03:30 committed by GitHub
parent 3d6fa28221 2ec3dfdc7e
commit 164eedb27c
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
27 changed files with 987 additions and 309 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -6,7 +6,7 @@ project(phasicFlow VERSION 1.0 )
set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_INSTALL_PREFIX ${phasicFlow_SOURCE_DIR} CACHE PATH "Install path of phasicFlow" FORCE)
set(CMAKE_BUILD_TYPE Debug CACHE STRING "build type" FORCE)
set(CMAKE_BUILD_TYPE Release CACHE STRING "build type" FORCE)
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build using shared libraries" FORCE)
mark_as_advanced(FORCE var BUILD_SHARED_LIBS)

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

@ -37,7 +37,8 @@ bool intAllActive(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1)
realx3PointField_D& dy1,
real damping = 1.0)
{
auto d_dy = dy.deviceView();
@ -49,7 +50,7 @@ bool intAllActive(
"AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] += dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[i]);
d_y[i] = damping*(d_dy[i] + 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();
@ -62,7 +63,8 @@ bool intScattered
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1
realx3PointField_D& dy1,
real damping = 1.0
)
{
@ -78,7 +80,7 @@ bool intScattered
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] += dt*(static_cast<real>(1.5) * d_dy[i] - static_cast<real>(0.5) * d_dy1[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_dy1[i] = d_dy[i];
}
});
@ -142,18 +144,19 @@ bool pFlow::AdamsBashforth2::correct
(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy
realx3PointField_D& dy,
real damping
)
{
auto& dy1l = dy1();
bool success = false;
if(dy1l.isAllActive())
{
success = intAllActive(dt, y.field(), dy, dy1l);
success = intAllActive(dt, y.field(), dy, dy1l, damping);
}
else
{
success = intScattered(dt, y.field(), dy, dy1l);
success = intScattered(dt, y.field(), dy, dy1l, damping);
}
success = success && boundaryList_.correct(dt, y, dy);

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

@ -45,6 +45,8 @@ private:
friend class processorAB2BoundaryIntegration;
protected:
const auto& dy1()const
{
return static_cast<const realx3PointField_D&>(*this);
@ -54,6 +56,11 @@ private:
{
return static_cast<realx3PointField_D&>(*this);
}
boundaryIntegrationList& boundaryList()
{
return boundaryList_;
}
public:
@ -70,7 +77,7 @@ public:
const realx3Field_D& initialValField);
/// Destructor
~AdamsBashforth2()final = default;
~AdamsBashforth2()override = default;
/// Add this to the virtual constructor table
add_vCtor(
@ -102,12 +109,13 @@ public:
bool correct(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy) final;
realx3PointField_D& dy,
real damping = 1.0) override;
bool correctPStruct(
real dt,
pointStructure& pStruct,
realx3PointField_D& vel) final;
realx3PointField_D& vel) override;
/*bool hearChanges
@ -121,9 +129,9 @@ public:
bool setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y) final;
const realx3Vector& y) override;
bool needSetInitialVals()const final
bool needSetInitialVals()const override
{
return false;
}

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

@ -20,62 +20,171 @@ Licence:
#include "AdamsBashforth3.hpp"
namespace pFlow
{
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool intAllActive(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
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 activeRng = dy1.activeRange();
Kokkos::parallel_for(
"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]
- static_cast<real>(16.0 / 12.0) * d_dy1[i]
+ static_cast<real>(5.0 / 12.0) * 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,
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 activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth2::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]
- static_cast<real>(16.0 / 12.0) * d_dy1[i]
+ static_cast<real>(5.0 / 12.0) * d_dy2[i]));
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
//const real AB3_coef[] = { 23.0 / 12.0, 16.0 / 12.0, 5.0 / 12.0 };
pFlow::AdamsBashforth3::AdamsBashforth3
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
integration(baseName, owner, pStruct, method),
history_(
owner.emplaceObject<pointField<VectorSingle,AB3History>>(
objectFile(
groupNames(baseName,"AB3History"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
AB3History({zero3,zero3})))
AdamsBashforth2(baseName, pStruct, method, initialValField),
dy2_
(
objectFile
(
groupNames(baseName,"dy2"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{
}
bool pFlow::AdamsBashforth3::predict
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
void pFlow::AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested()
{
return true;
AdamsBashforth2::updateBoundariesSlaveToMasterIfRequested();
dy2_.updateBoundariesSlaveToMasterIfRequested();
}
bool pFlow::AdamsBashforth3::correct
(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
real dt,
realx3PointField_D& y,
realx3PointField_D& dy,
real damping
)
{
if(this->pStruct().allActive())
bool success = false;
if(y.isAllActive())
{
return intAll(dt, y, dy, this->pStruct().activeRange());
success = intAllActive(dt, y.field(), dy, dy1(), dy2(), damping);
}
else
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
success = intScattered(dt, y.field(), dy, dy1(), dy2(), damping);
}
return true;
success = success && boundaryList().correct(dt, y, dy);
return success;
}
bool pFlow::AdamsBashforth3::correctPStruct(
real dt,
pointStructure &pStruct,
realx3PointField_D &vel)
{
bool success = false;
if(dy2().isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2());
}
success = success && boundaryList().correctPStruct(dt, pStruct, vel);
return success;
}
bool pFlow::AdamsBashforth3::setInitialVals(
const int32IndexContainer& newIndices,
const realx3Vector& y)
@ -83,7 +192,7 @@ bool pFlow::AdamsBashforth3::setInitialVals(
return true;
}
bool pFlow::AdamsBashforth3::intAll(
/*bool pFlow::AdamsBashforth3::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
@ -106,4 +215,4 @@ bool pFlow::AdamsBashforth3::intAll(
Kokkos::fence();
return true;
}
}*/

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

@ -22,48 +22,14 @@ Licence:
#define __AdamsBashforth3_hpp__
#include "integration.hpp"
#include "AdamsBashforth2.hpp"
#include "pointFields.hpp"
#include "boundaryIntegrationList.hpp"
namespace pFlow
{
struct AB3History
{
TypeInfoNV("AB3History");
realx3 dy1_={0,0,0};
realx3 dy2_={0,0,0};
};
INLINE_FUNCTION
iIstream& operator>>(iIstream& str, AB3History& ab3)
{
str.readBegin("AB3History");
str >> ab3.dy1_;
str >> ab3.dy2_;
str.readEnd("AB3History");
str.check(FUNCTION_NAME);
return str;
}
INLINE_FUNCTION
iOstream& operator<<(iOstream& str, const AB3History& ab3)
{
str << token::BEGIN_LIST << ab3.dy1_
<< token::SPACE << ab3.dy2_
<< token::END_LIST;
str.check(FUNCTION_NAME);
return str;
}
/**
* Third order Adams-Bashforth integration method for solving ODE
@ -72,19 +38,26 @@ iOstream& operator<<(iOstream& str, const AB3History& ab3)
*/
class AdamsBashforth3
:
public integration
public AdamsBashforth2
{
protected:
private:
/// Integration history
pointField<VectorSingle,AB3History>& history_;
friend class processorAB3BoundaryIntegration;
/// Range policy for integration kernel
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
realx3PointField_D dy2_;
protected:
const auto& dy2()const
{
return dy2_;
}
auto& dy2()
{
return dy2_;
}
public:
@ -96,17 +69,13 @@ public:
/// Construct from components
AdamsBashforth3(
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<AdamsBashforth3>(*this);
}
/// Destructor
virtual ~AdamsBashforth3()=default;
~AdamsBashforth3() override =default;
/// Add this to the virtual constructor table
add_vCtor(
@ -117,14 +86,36 @@ public:
// - Methods
bool predict(
real UNUSED(dt),
realx3Vector_D & UNUSED(y),
realx3Vector_D& UNUSED(dy)) override;
void updateBoundariesSlaveToMasterIfRequested()override;
bool correct(real dt,
realx3Vector_D & y,
realx3Vector_D& dy) override;
/// return integration method
word method()const override
{
return "AdamsBashforth3";
}
bool correct(
real dt,
realx3PointField_D& y,
realx3PointField_D& dy,
real damping = 1.0) override;
bool correctPStruct(
real dt,
pointStructure& pStruct,
realx3PointField_D& vel) override;
/*bool hearChanges
(
real t,
real dt,
uint32 iter,
const message& msg,
const anyList& varList
) override;*/
bool setInitialVals(
const int32IndexContainer& newIndices,
@ -135,25 +126,10 @@ public:
return false;
}
/// Integrate on all points in the active range
bool intAll(
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 );
};
template<typename activeFunctor>
/*template<typename activeFunctor>
bool pFlow::AdamsBashforth3::intRange(
real dt,
realx3Vector_D& y,
@ -181,7 +157,7 @@ bool pFlow::AdamsBashforth3::intRange(
Kokkos::fence();
return true;
}
}*/
} // pFlow

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

@ -20,60 +20,174 @@ Licence:
#include "AdamsBashforth4.hpp"
namespace pFlow
{
/// Range policy for integration kernel (alias)
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<uint32>
>;
bool intAllActive(
real dt,
realx3Field_D& y,
realx3PointField_D& dy,
realx3PointField_D& dy1,
realx3PointField_D& dy2,
realx3PointField_D& dy3,
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 activeRng = dy1.activeRange();
Kokkos::parallel_for(
"AdamsBashforth3::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_y[i] = damping *(d_y[i] + 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]
));
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,
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 activeRng = dy1.activeRange();
const auto& activeP = dy1.activePointsMaskDevice();
Kokkos::parallel_for(
"AdamsBashforth2::correct",
rpIntegration (activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
if( activeP(i))
{
d_y[i] = damping* ( d_y[i] + 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]
));
d_dy3[i] = d_dy2[i];
d_dy2[i] = d_dy1[i];
d_dy1[i] = d_dy[i];
}
});
Kokkos::fence();
return true;
}
}
pFlow::AdamsBashforth4::AdamsBashforth4
(
const word& baseName,
repository& owner,
const pointStructure& pStruct,
const word& method
pointStructure& pStruct,
const word& method,
const realx3Field_D& initialValField
)
:
integration(baseName, owner, pStruct, method),
history_(
owner.emplaceObject<pointField<VectorSingle,AB4History>>(
objectFile(
groupNames(baseName,"AB4History"),
"",
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS),
pStruct,
AB4History({zero3,zero3, zero3})))
AdamsBashforth3(baseName, pStruct, method, initialValField),
dy3_
(
objectFile
(
groupNames(baseName,"dy3"),
pStruct.time().integrationFolder(),
objectFile::READ_IF_PRESENT,
objectFile::WRITE_ALWAYS
),
pStruct,
zero3,
zero3
)
{
}
bool pFlow::AdamsBashforth4::predict
(
real UNUSED(dt),
realx3Vector_D& UNUSED(y),
realx3Vector_D& UNUSED(dy)
)
void pFlow::AdamsBashforth4::updateBoundariesSlaveToMasterIfRequested()
{
return true;
AdamsBashforth3::updateBoundariesSlaveToMasterIfRequested();
dy3_.updateBoundariesSlaveToMasterIfRequested();
}
bool pFlow::AdamsBashforth4::correct
(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy
real dt,
realx3PointField_D& y,
realx3PointField_D& dy,
real damping
)
{
if(this->pStruct().allActive())
bool success = false;
if(y.isAllActive())
{
return intAll(dt, y, dy, this->pStruct().activeRange());
success = intAllActive(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
}
else
{
return intRange(dt, y, dy, this->pStruct().activePointsMaskD());
success = intScattered(dt, y.field(), dy, dy1(), dy2(), dy3(), damping);
}
return true;
success = success && boundaryList().correct(dt, y, dy);
return success;
}
bool pFlow::AdamsBashforth4::correctPStruct(
real dt,
pointStructure &pStruct,
realx3PointField_D &vel)
{
bool success = false;
if(dy2().isAllActive())
{
success = intAllActive(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
}
else
{
success = intScattered(dt, pStruct.pointPosition(), vel, dy1(), dy2(), dy3());
}
success = success && boundaryList().correctPStruct(dt, pStruct, vel);
return success;
}
bool pFlow::AdamsBashforth4::setInitialVals(
@ -83,7 +197,7 @@ bool pFlow::AdamsBashforth4::setInitialVals(
return true;
}
bool pFlow::AdamsBashforth4::intAll(
/*bool pFlow::AdamsBashforth4::intAll(
real dt,
realx3Vector_D& y,
realx3Vector_D& dy,
@ -112,4 +226,4 @@ bool pFlow::AdamsBashforth4::intAll(
Kokkos::fence();
return true;
}
}*/

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

@ -22,53 +22,14 @@ Licence:
#define __AdamsBashforth4_hpp__
#include "integration.hpp"
#include "AdamsBashforth3.hpp"
#include "pointFields.hpp"
#include "boundaryIntegrationList.hpp"
namespace pFlow
{
struct AB4History
{
TypeInfoNV("AB4History");
realx3 dy1_={0,0,0};
realx3 dy2_={0,0,0};
realx3 dy3_={0,0,0};
};
INLINE_FUNCTION
iIstream& operator>>(iIstream& str, AB4History& ab4)
{
str.readBegin("AB4History");
str >> ab4.dy1_;
str >> ab4.dy2_;
str >> ab4.dy3_;
str.readEnd("AB4History");
str.check(FUNCTION_NAME);
return str;
}
INLINE_FUNCTION
iOstream& operator<<(iOstream& str, const AB4History& ab4)
{
str << token::BEGIN_LIST << ab4.dy1_
<< token::SPACE << ab4.dy2_
<< token::SPACE << ab4.dy3_
<< token::END_LIST;
str.check(FUNCTION_NAME);
return str;
}
/**
* Fourth order Adams-Bashforth integration method for solving ODE
*
@ -76,19 +37,25 @@ iOstream& operator<<(iOstream& str, const AB4History& ab4)
*/
class AdamsBashforth4
:
public integration
public AdamsBashforth3
{
private:
friend class processorAB4BoundaryIntegration;
realx3PointField_D dy3_;
protected:
/// Integration history
pointField<VectorSingle,AB4History>& history_;
const auto& dy3()const
{
return dy3_;
}
/// Range policy for integration kernel
using rpIntegration = Kokkos::RangePolicy<
DefaultExecutionSpace,
Kokkos::Schedule<Kokkos::Static>,
Kokkos::IndexType<int32>
>;
auto& dy3()
{
return dy3_;
}
public:
@ -100,17 +67,14 @@ public:
/// Construct from components
AdamsBashforth4(
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<AdamsBashforth4>(*this);
}
/// Destructor
virtual ~AdamsBashforth4()=default;
~AdamsBashforth4() override =default;
/// Add a this to the virtual constructor table
add_vCtor(
@ -121,15 +85,24 @@ 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 "AdamsBashforth4";
}
bool correct(
real dt,
realx3Vector_D & y,
realx3Vector_D& dy) override;
realx3PointField_D& y,
realx3PointField_D& dy,
real damping = 1.0) override;
bool correctPStruct(
real dt,
pointStructure& pStruct,
realx3PointField_D& vel) override;
bool setInitialVals(
const int32IndexContainer& newIndices,
@ -140,25 +113,12 @@ public:
return false;
}
/// Integrate on all points in the active range
bool intAll(
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 );
};
template<typename activeFunctor>
/*template<typename activeFunctor>
bool pFlow::AdamsBashforth4::intRange(
real dt,
realx3Vector_D& y,
@ -191,7 +151,7 @@ bool pFlow::AdamsBashforth4::intRange(
Kokkos::fence();
return true;
}
}*/
} // pFlow

5
src/Integration/CMakeLists.txt
View File

@ -4,9 +4,10 @@ integration/integration.cpp
boundaries/boundaryIntegration.cpp
boundaries/boundaryIntegrationList.cpp
AdamsBashforth2/AdamsBashforth2.cpp
AdamsBashforth3/AdamsBashforth3.cpp
AdamsBashforth4/AdamsBashforth4.cpp
#AdamsBashforth5/AdamsBashforth5.cpp
#AdamsBashforth4/AdamsBashforth4.cpp
#AdamsBashforth3/AdamsBashforth3.cpp
#AdamsMoulton3/AdamsMoulton3.cpp
#AdamsMoulton4/AdamsMoulton4.cpp
#AdamsMoulton5/AdamsMoulton5.cpp

2
src/Integration/integration/integration.hpp
View File

@ -146,7 +146,7 @@ public:
/// Correction/main integration step
virtual
bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy) = 0;
bool correct(real dt, realx3PointField_D& y, realx3PointField_D& dy, real damping = 1.0) = 0;
virtual
bool correctPStruct(real dt, pointStructure& pStruct, realx3PointField_D& vel) = 0;

307
src/Interaction/Models/contactForce/cGNonLinearCF.hpp Normal file
View File

@ -0,0 +1,307 @@
/*------------------------------- 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 __cGNonLinearCF_hpp__
#define __cGNonLinearCF_hpp__
#include "types.hpp"
namespace pFlow::cfModels
{
template<bool limited=true>
class cGNonLinear
{
public:
struct contactForceStorage
{
realx3 overlap_t_ = 0.0;
};
struct cGNonLinearProperties
{
real Yeff_ = 1000000.0;
real Geff_ = 8000000.0;
real en_ = 1.0;
real mu_ = 0.00001;
INLINE_FUNCTION_HD
cGNonLinearProperties(){}
INLINE_FUNCTION_HD
cGNonLinearProperties(real Yeff, real Geff, real en, real mu ):
Yeff_(Yeff), Geff_(Geff), en_(en), mu_(mu)
{}
INLINE_FUNCTION_HD
cGNonLinearProperties(const cGNonLinearProperties&)=default;
INLINE_FUNCTION_HD
cGNonLinearProperties& operator=(const cGNonLinearProperties&)=default;
INLINE_FUNCTION_HD
~cGNonLinearProperties() = default;
};
protected:
using cGNonLinearArrayType = symArray<cGNonLinearProperties>;
int32 numMaterial_ = 0;
ViewType1D<real> rho_;
int32 addDissipationModel_;
cGNonLinearArrayType nonlinearProperties_;
bool readNonLinearDictionary(const dictionary& dict)
{
auto Yeff = dict.getVal<realVector>("Yeff");
auto Geff = dict.getVal<realVector>("Geff");
auto nu = dict.getVal<realVector>("nu");
auto en = dict.getVal<realVector>("en");
auto mu = dict.getVal<realVector>("mu");
auto nElem = Yeff.size();
if(nElem != nu.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and nu("<<nu.size()<<") do not match.\n";
return false;
}
if(nElem != en.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and en("<<en.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size())
{
fatalErrorInFunction<<
"sizes of Yeff("<<nElem<<") and mu("<<mu.size()<<") do not match.\n";
return false;
}
// check if size of vector matchs a symetric array
uint32 nMat;
if( !cGNonLinearArrayType::getN(nElem, nMat) )
{
fatalErrorInFunction<<
"sizes of properties do not match a symetric array with size ("<<
numMaterial_<<"x"<<numMaterial_<<").\n";
return false;
}
else if( numMaterial_ != nMat)
{
fatalErrorInFunction<<
"size mismatch for porperties. \n";
return false;
}
Vector<cGNonLinearProperties> prop("prop",nElem);
ForAll(i,Yeff)
{
prop[i] = {Yeff[i], Geff[i], en[i], mu[i]};
}
nonlinearProperties_.assign(prop);
auto adm = dict.getVal<word>("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 "cGNonLinearLimited";
}
else
{
return "cGNonLinearNonLimited";
}
return "";
}
public:
TypeInfoNV(modelName());
INLINE_FUNCTION_HD
cGNonLinear(){}
cGNonLinear(
int32 nMaterial,
const ViewType1D<real>& rho,
const dictionary& dict)
:
numMaterial_(nMaterial),
rho_("rho",nMaterial),
nonlinearProperties_("cGNonLinearProperties",nMaterial)
{
Kokkos::deep_copy(rho_,rho);
if(!readNonLinearDictionary(dict))
{
fatalExit;
}
}
INLINE_FUNCTION_HD
cGNonLinear(const cGNonLinear&) = default;
INLINE_FUNCTION_HD
cGNonLinear(cGNonLinear&&) = default;
INLINE_FUNCTION_HD
cGNonLinear& operator=(const cGNonLinear&) = default;
INLINE_FUNCTION_HD
cGNonLinear& operator=(cGNonLinear&&) = default;
INLINE_FUNCTION_HD
~cGNonLinear()=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
{
const auto prop = nonlinearProperties_(propId_i,propId_j);
const real f = 2.0/( 1.0/cGFi + 1.0/cGFj );
real vrn = dot(Vr, Nij);
realx3 Vt = Vr - vrn*Nij;
history.overlap_t_ += Vt*dt;
real mi = 3*Pi/4*pow(Ri,static_cast<real>(3))*rho_[propId_i];
real mj = 3*Pi/4*pow(Rj,static_cast<real>(3))*rho_[propId_j];
real Reff = 1.0/(1/Ri + 1/Rj);
real K_hertz = 4.0/3.0*prop.Yeff_*sqrt(Reff);
real sqrt_meff_K_hertz = sqrt((mi*mj)/(mi+mj) * K_hertz);
real en = prop.en_;
if (addDissipationModel_==2)
{
en = sqrt(1+((pow(prop.en_,2)-1)*f));
}
else if (addDissipationModel_==3)
{
en = exp((pow(f,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
}
real Kn = static_cast<real>(4.0/3.0) * prop.Yeff_ * sqrt(Reff*ovrlp_n);
real ethan_ = -2.0*log(en)*sqrt(Kn)/
sqrt(pow(log(en),2.0)+ pow(Pi,2.0));
FCn = ( - Kn*ovrlp_n -
sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast<real>(0.25))*vrn)*Nij;
FCt = (f*f)*(- static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n) ) * history.overlap_t_;
real ft = length(FCt);
real ft_fric = prop.mu_ * length(FCn);
// apply friction
if(ft > ft_fric)
{
if( length(history.overlap_t_) >0.0)
{
if constexpr (limited)
{
real kt = static_cast<real>(8.0) * prop.Geff_ * sqrt(Reff*ovrlp_n);
FCt *= (ft_fric/ft);
history.overlap_t_ = - FCt/(f*f*kt);
}
else
{
FCt = (FCt/ft)*ft_fric;
}
}
else
{
FCt = 0.0;
}
}
}
};
} //pFlow::CFModels
#endif

60
src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp
View File

@ -26,11 +26,6 @@ Licence:
namespace pFlow::cfModels
{
@ -49,16 +44,15 @@ public:
{
real kn_ = 1000.0;
real kt_ = 800.0;
real en_ = 0.0;
real ethat_ = 0.0;
real mu_ = 0.00001;
real en_ = 1.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)
linearProperties(real kn, real kt, real en, real mu ):
kn_(kn), kt_(kt), en_(en), mu_(mu)
{}
INLINE_FUNCTION_HD
@ -93,7 +87,6 @@ protected:
auto kn = dict.getVal<realVector>("kn");
auto kt = dict.getVal<realVector>("kt");
auto en = dict.getVal<realVector>("en");
auto et = dict.getVal<realVector>("et");
auto mu = dict.getVal<realVector>("mu");
@ -114,12 +107,6 @@ protected:
return false;
}
if(nElem != et.size())
{
fatalErrorInFunction<<
"sizes of kn("<<nElem<<") and et("<<et.size()<<") do not match.\n";
return false;
}
if(nElem != mu.size())
{
@ -148,30 +135,16 @@ protected:
}
realVector etha_t("etha_t", nElem);
ForAll(i , kn)
{
etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) /
sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0));
}
Vector<linearProperties> prop("prop", nElem);
ForAll(i,kn)
{
prop[i] = {kn[i], kt[i], en[i], etha_t[i], mu[i] };
prop[i] = {kn[i], kt[i], en[i], mu[i] };
}
linearProperties_.assign(prop);
auto adm = dict.getVal<word>("additionalDissipationModel");
auto adm = dict.getVal<word>("additionalDissipationModel");
if(adm == "none")
{
addDissipationModel_ = 1;
@ -287,22 +260,23 @@ public:
real sqrt_meff = sqrt((mi*mj)/(mi+mj));
if (addDissipationModel_==2)
real en = prop.en_;
if (addDissipationModel_==2)
{
prop.en_ = sqrt(1+((pow(prop.en_,2)-1)*f_));
en = sqrt(1+((pow(prop.en_,2)-1)*f_));
}
else if (addDissipationModel_==3)
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)))) ) ));
en = exp((pow(f_,1.5)*log(prop.en_)*sqrt( (1-((pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2))))/(1-(pow(f_,3)*(pow(log(prop.en_),2))/(pow(log(prop.en_),2)+pow(Pi,2)))) ) ));
}
real ethan_ = -2.0*log(prop.en_)*sqrt(prop.kn_)/
sqrt(pow(log(prop.en_),2.0)+ pow(Pi,2.0));
real ethan_ = -2.0*log(en)*sqrt(prop.kn_)/
sqrt(pow(log(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);
FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij;
FCt = ( -f_*prop.kt_ * history.overlap_t_);

4
src/Interaction/Models/grainContactForceModels.hpp
View File

@ -23,6 +23,7 @@ Licence:
#include "cGAbsoluteLinearCF.hpp"
#include "cGRelativeLinearCF.hpp"
#include "cGNonLinearCF.hpp"
#include "grainRolling.hpp"
@ -38,6 +39,9 @@ using nonLimitedCGAbsoluteLinearGrainRolling = grainRolling<cGAbsoluteLinear<fal
using limitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<true>>;
using nonLimitedCGRelativeLinearGrainRolling = grainRolling<cGRelativeLinear<false>>;
using limitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<true>>;
using nonLimitedCGNonLinearGrainRolling = grainRolling<cGNonLinear<false>>;
}

6
src/Interaction/grainInteraction/grainInteractionsNonLinearModels.cpp
View File

@ -41,3 +41,9 @@ Licence:
createBoundaryGrainInteraction(ForceModel, GeomModel)
createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::rotationAxisMotionGeometry);
createInteraction(pFlow::cfModels::limitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);
createInteraction(pFlow::cfModels::nonLimitedCGNonLinearGrainRolling, pFlow::stationaryGeometry);

3
src/Particles/CMakeLists.txt
View File

@ -6,6 +6,9 @@ particles/shape/shape.cpp
particles/particles.cpp
particles/particleIdHandler/particleIdHandler.cpp
particles/regularParticleIdHandler/regularParticleIdHandler.cpp
globalDamping/globalDamping.cpp
SphereParticles/sphereShape/sphereShape.cpp
SphereParticles/sphereParticles/sphereParticles.cpp
SphereParticles/sphereParticles/sphereParticlesKernels.cpp

17
src/Particles/dynamicPointStructure/dynamicPointStructure.cpp
View File

@ -21,6 +21,7 @@ Licence:
#include "dynamicPointStructure.hpp"
#include "systemControl.hpp"
pFlow::dynamicPointStructure::dynamicPointStructure
(
systemControl& control
@ -77,6 +78,9 @@ pFlow::dynamicPointStructure::dynamicPointStructure
fatalExit;
}
REPORT(1)<<"Reading globalDamping dictionary ..."<<END_REPORT;
velDamping_ = makeUnique<globalDamping>(control);
}
@ -101,6 +105,16 @@ bool pFlow::dynamicPointStructure::iterate()
return correct(dt, acc);*/
}
bool pFlow::dynamicPointStructure::afterIteration()
{
//const auto ti = TimeInfo();
auto succs = pointStructure::afterIteration();
//velDamping_().applyDamping(ti, velocity_);
return succs;
}
bool pFlow::dynamicPointStructure::predict(
real dt,
realx3PointField_D &acceleration)
@ -119,10 +133,11 @@ bool pFlow::dynamicPointStructure::correct
)
{
//auto& pos = pStruct().pointPosition();
const auto ti = TimeInfo();
if(!integrationPos_().correctPStruct(dt, *this, velocity_) )return false;
if(!integrationVel_().correct(dt, velocity_, acceleration))return false;
if(!integrationVel_().correct(dt, velocity_, acceleration, velDamping_().dampingFactor(ti)))return false;
return true;
}

6
src/Particles/dynamicPointStructure/dynamicPointStructure.hpp
View File

@ -26,11 +26,13 @@ Licence:
#include "pointFields.hpp"
#include "integration.hpp"
#include "uniquePtr.hpp"
#include "globalDamping.hpp"
namespace pFlow
{
class systemControl;
//class globalDamping;
class dynamicPointStructure
:
@ -44,6 +46,8 @@ private:
uniquePtr<integration> integrationVel_ = nullptr;
uniquePtr<globalDamping> velDamping_ = nullptr;
Timer velocityUpdateTimer_;
/// @brief integration method for velocity and position
@ -88,6 +92,8 @@ public:
/// @brief This is called in time loop. Perform the main calculations
/// when the component should evolve along time.
bool iterate() override;
bool afterIteration()override;
/// prediction step (if any), is called in beforeIteration
bool predict(real dt, realx3PointField_D& acceleration);

74
src/Particles/globalDamping/globalDamping.cpp Normal file
View File

@ -0,0 +1,74 @@
/*------------------------------- phasicFlow ---------------------------------
O C enter of
O O E ngineering and
O O M ultiscale modeling of
OOOOOOO F luid flow
------------------------------------------------------------------------------
Copyright (C): www.cemf.ir
email: hamid.r.norouzi AT gmail.com
------------------------------------------------------------------------------
Licence:
This file is part of phasicFlow code. It is a free software for simulating
granular and multiphase flows. You can redistribute it and/or modify it under
the terms of GNU General Public License v3 or any other later versions.
phasicFlow is distributed to help others in their research in the field of
granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the
implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-----------------------------------------------------------------------------*/
#include "globalDamping.hpp"
pFlow::globalDamping::globalDamping(const systemControl& control)
:
timeControl_(control.settingsDict().subDict("globalDamping"), control.time().dt(), "damping")
{
const dictionary& dict = control.settingsDict().subDict("globalDamping");
dampingFactor_ = dict.getValOrSetMin<real>("dampingFactor", static_cast<real>(1.0));
dampingFactor_ = max( dampingFactor_ , static_cast<real>(0.01));
performDamping_ = !equal(dampingFactor_, static_cast<real>(1.0));
if( performDamping_ )
REPORT(2)<<"Global damping "<<Yellow_Text("is active")<<
" and damping factor is "<<dampingFactor_<<END_REPORT;
else
REPORT(2)<<"Global damping "<<Yellow_Text("is not active")<<"."<<END_REPORT;
}
/*void pFlow::globalDamping::applyDamping
(
const timeInfo& ti,
realx3PointField_D& velocity
)
{
if(!performDamping_) return;
if(!timeControl_.timeEvent(ti.iter(), ti.t(), ti.dt()) )return;
auto d_v = velocity.deviceView();
auto activeRng = velocity.activeRange();
auto dmpng = dampingFactor_;
Kokkos::parallel_for(
"globalDamping::applyDamping",
deviceRPolicyStatic(activeRng.start(), activeRng.end()),
LAMBDA_HD(uint32 i){
d_v[i] *= dmpng;
});
Kokkos::fence();
//REPORT(1)<<"Applied global damping "<<END_REPORT;
}*/
pFlow::real pFlow::globalDamping::dampingFactor(const timeInfo& ti)const
{
if(!performDamping_) return 1.0;
if(!timeControl_.timeEvent(ti.iter(), ti.t(), ti.dt()) )return 1.0;
return dampingFactor_;
}

65
src/Particles/globalDamping/globalDamping.hpp Normal file
View File

@ -0,0 +1,65 @@
/*------------------------------- 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 __globalDamping_hpp__
#define __globalDamping_hpp__
#include "systemControl.hpp"
#include "pointFields.hpp"
#include "baseTimeControl.hpp"
namespace pFlow
{
class globalDamping
{
private:
bool performDamping_ = false;
real dampingFactor_;
baseTimeControl timeControl_;
public:
globalDamping(const systemControl& control);
~globalDamping()=default;
//void applyDamping( const timeInfo& ti, realx3PointField_D& velocity);
bool performDamping()const
{
return performDamping_;
}
real dampingFactor(const timeInfo& ti)const;
};
}
#endif

2
src/phasicFlow/dictionary/dictionary.hpp
View File

@ -431,7 +431,7 @@ T dictionary::getValOrSet
template<typename T>
T dictionary::getValOrSetMax(const word& keyword, const T& setMaxVal)const
{
return max(getValOrSet(keyword, setMaxVal), setMaxVal);
return max(getValOrSet(keyword, setMaxVal), static_cast<T>(setMaxVal));
}
template<typename T>

44
src/phasicFlow/repository/Time/baseTimeControl.cpp
View File

@ -65,6 +65,50 @@ pFlow::baseTimeControl::baseTimeControl
}
pFlow::baseTimeControl::baseTimeControl
(
const dictionary& dict,
const real defInterval,
const word& intervalPrefix,
const real defStartTime
)
:
intervalPrefix_(intervalPrefix)
{
auto tControl = dict.getVal<word>("timeControl");
if(tControl == "timeStep")
isTimeStep_ = true;
else if( tControl == "runTime" ||
tControl == "simulationTime")
isTimeStep_ = false;
else
{
fatalErrorInFunction<<
"bad input for intervalControl in dictionary "<< dict.globalName()<<endl<<
"valid inputs are "<<
wordList({"timeStep", "runTime", "simulationTime"})<<endl;
fatalExit;
}
word intervalWord = intervalPrefix.size()==0? word("interval"): intervalPrefix+"Interval";
if(!isTimeStep_)
{
auto startTime = (dict.getValOrSet<real>("startTime", defStartTime));
auto endTime = (dict.getValOrSet<real>("endTime", largeValue));
auto interval = dict.getValOrSet<real>(intervalWord, defInterval);
rRange_ = realStridedRange(startTime, endTime, interval);
}
else
{
auto startTime = (dict.getValOrSet<int32>("startTime", 0));
auto endTime = (dict.getValOrSet<int32>("endTime", largestPosInt32));
auto interval = dict.getValOrSet<int32>(intervalWord, 1);
iRange_ = int32StridedRagne(startTime, endTime, interval);
}
}
pFlow::baseTimeControl::baseTimeControl(int32 start, int32 end, int32 stride, const word &intervalPrefix)
:
isTimeStep_(true),

6
src/phasicFlow/repository/Time/baseTimeControl.hpp
View File

@ -46,6 +46,12 @@ public:
real defStartTime = 0.0
);
baseTimeControl(
const dictionary& dict,
const real defInterval,
const word& intervalPrefix="",
const real defStartTime=0.0);
baseTimeControl(
int32 start,
int32 end,

2
src/phasicFlow/repository/Time/timeControl.cpp
View File

@ -167,7 +167,7 @@ void pFlow::timeControl::checkForOutputToFile()
lastSaved_ = currentTime_;
save = true;
}
else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), 0.5 *dt_) )
else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), static_cast<real>(0.5 *dt_)) )
{
lastSaved_ = currentTime_;
save = true;

2
src/phasicFlow/structuredData/domain/simulationDomain.cpp
View File

@ -31,7 +31,7 @@ bool pFlow::simulationDomain::prepareBoundaryDicts()
real neighborLength = boundaries.getVal<real>("neighborLength");
real boundaryExtntionLengthRatio =
boundaries.getValOrSetMax("boundaryExtntionLengthRatio", 0.1);
boundaries.getValOrSetMax("boundaryExtntionLengthRatio", static_cast<real>(0.1));
uint32 updateInterval =
boundaries.getValOrSetMax<uint32>("updateInterval", 1u);

2
src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp
View File

@ -28,7 +28,7 @@ pFlow::infinitePlane::infinitePlane
)
{
auto ln = cross(p2-p1, p3-p1);
if( equal(ln.length(),0.0) )
{
fatalErrorInFunction<<

1
src/phasicFlow/types/basicTypes/bTypesFunctions.cpp
View File

@ -21,7 +21,6 @@ Licence:
#include <algorithm>
#include <iomanip>
#include <sstream>
#include "bTypesFunctions.hpp"
pFlow::int32

2
src/phasicFlow/types/basicTypes/bTypesFunctions.hpp
View File

@ -26,7 +26,7 @@ Licence:
#define __bTypesFunctions_hpp__
#include "builtinTypes.hpp"
//#include "math.hpp"
#include "math.hpp"
#include "numericConstants.hpp"
#include "pFlowMacros.hpp"

50
utilities/particlesPhasicFlow/particlesPhasicFlow.cpp
View File

@ -181,32 +181,36 @@ int main( int argc, char* argv[] )
&Control.caseSetup()
);
auto& shapeName = Control.time().template lookupObject<pFlow::wordPointField_H>("shapeName");
REPORT(0)<< "Converting shapeName field to shapeIndex field"<<END_REPORT;
auto shapeName_D = shapeName.deviceView();
auto shapeHash_D = shapeHash.deviceView();
auto shapeIndex_D = shapeIndex.deviceView();
REPORT(1)<<"List of shape names in "<<shapes.globalName()<<
" is: "<<Green_Text(shapes.shapeNameList())<<END_REPORT;
ForAll(i, shapeHash)
if(Control.time().lookupObjectName("shapeName"))
{
if(pFlow::uint32 index; shapes.shapeNameToIndex(shapeName_D[i], index))
auto& shapeName = Control.time().template lookupObject<pFlow::wordPointField_H>("shapeName");
REPORT(0)<< "Converting shapeName field to shapeIndex field"<<END_REPORT;
auto shapeName_D = shapeName.deviceView();
auto shapeHash_D = shapeHash.deviceView();
auto shapeIndex_D = shapeIndex.deviceView();
REPORT(1)<<"List of shape names in "<<shapes.globalName()<<
" is: "<<Green_Text(shapes.shapeNameList())<<END_REPORT;
ForAll(i, shapeHash)
{
shapeHash_D[i] = shapes.hashes()[index];
shapeIndex_D[i] = index;
}
else
{
fatalErrorInFunction<<"Found shape name "<< Yellow_Text(shapeName_D[i])<<
"in shapeName field. But the list of shape names in file "<<
shapes.globalName()<<" is : \n"<<
shapes.shapeNames()<<pFlow::endl;
}
if(pFlow::uint32 index; shapes.shapeNameToIndex(shapeName_D[i], index))
{
shapeHash_D[i] = shapes.hashes()[index];
shapeIndex_D[i] = index;
}
else
{
fatalErrorInFunction<<"Found shape name "<< Yellow_Text(shapeName_D[i])<<
"in shapeName field. But the list of shape names in file "<<
shapes.globalName()<<" is : \n"<<
shapes.shapeNames()<<pFlow::endl;
}
}
}
if( !Control.time().write(true))
{