diff --git a/src/Integration/AdamsMoulton3/AdamsMoulton3.C b/src/Integration/AdamsMoulton3/AdamsMoulton3.C index aa6d45d9..7c7cbe58 100644 --- a/src/Integration/AdamsMoulton3/AdamsMoulton3.C +++ b/src/Integration/AdamsMoulton3/AdamsMoulton3.C @@ -113,7 +113,7 @@ bool pFlow::AdamsMoulton3::setInitialVals( const realx3Vector& y) { y0_.insertSetElement(newIndices, y); - //output<< "y0_ "<< y0_<(5.0/12.0)*d_dy[i] + static_cast(8.0/12.0)*d_dy0[i] - static_cast(1.0/12.0)*d_dy1[i]); - d_y[i] = corrct_y; - d_y0[i] = corrct_y; d_dy1[i]= d_dy0[i]; + d_y0[i] = corrct_y; + d_y[i] = corrct_y; } }); Kokkos::fence(); diff --git a/src/Integration/AdamsMoulton4/AdamsMoulton4.C b/src/Integration/AdamsMoulton4/AdamsMoulton4.C new file mode 100644 index 00000000..be3242e7 --- /dev/null +++ b/src/Integration/AdamsMoulton4/AdamsMoulton4.C @@ -0,0 +1,194 @@ +/*------------------------------- 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 "AdamsMoulton4.H" + +//const real AB2_coef[] = { 3.0 / 2.0, 1.0 / 2.0}; + +pFlow::AdamsMoulton4::AdamsMoulton4 +( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method +) +: + integration(baseName, owner, pStruct, method), + y0_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"y0"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3, + false + ) + ), + dy0_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"dy0"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3 + ) + ), + dy1_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"dy1"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3 + ) + ), + dy2_( + owner.emplaceObject( + objectFile( + groupNames(baseName,"dy2"), + "", + objectFile::READ_IF_PRESENT, + objectFile::WRITE_ALWAYS), + pStruct, + zero3 + ) + ) +{ + +} + +bool pFlow::AdamsMoulton4::predict +( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy +) +{ + + if(this->pStruct().allActive()) + { + return predictAll(dt, y, dy, this->pStruct().activeRange()); + } + else + { + return predictRange(dt, y, dy, this->pStruct().activePointsMaskD()); + } + + return true; +} + +bool pFlow::AdamsMoulton4::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::AdamsMoulton4::setInitialVals( + const int32IndexContainer& newIndices, + const realx3Vector& y) +{ + y0_.insertSetElement(newIndices, y); + + return true; +} + +bool pFlow::AdamsMoulton4::predictAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng) +{ + + auto d_dy = dy.deviceVectorAll(); + auto d_y = y.deviceVectorAll(); + + auto d_y0 = y0_.deviceVectorAll(); + auto d_dy0 = dy0_.deviceVectorAll(); + auto d_dy1 = dy1_.deviceVectorAll(); + auto d_dy2 = dy2_.deviceVectorAll(); + + Kokkos::parallel_for( + "AdamsMoulton4::predict", + rpIntegration (activeRng.first, activeRng.second), + LAMBDA_HD(int32 i){ + d_dy0[i] = d_dy[i]; + d_y[i] = d_y0[i] + dt*( + static_cast(23.0 /12.0 ) * d_dy[i] + - static_cast(16.0 / 12.0) * d_dy1[i] + + static_cast( 5.0 / 12.0) * d_dy2[i]); + }); + Kokkos::fence(); + + return true; +} + +bool pFlow::AdamsMoulton4::intAll( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + range activeRng) +{ + + auto d_dy = dy.deviceVectorAll(); + auto d_y = y.deviceVectorAll(); + + auto d_dy0 = dy0_.deviceVectorAll(); + auto d_y0 = y0_.deviceVectorAll(); + auto d_dy1 = dy1_.deviceVectorAll(); + auto d_dy2 = dy2_.deviceVectorAll(); + + Kokkos::parallel_for( + "AdamsMoulton4::correct", + rpIntegration (activeRng.first, activeRng.second), + LAMBDA_HD(int32 i){ + auto corrct_y = d_y0[i] + dt*( + static_cast(9.0/24.0)*d_dy[i] + + static_cast(19.0/24.0)*d_dy0[i] + - static_cast( 5.0/24.0)*d_dy1[i] + + static_cast( 1.0/24.0)*d_dy2[i]); + + d_dy2[i]= d_dy1[i]; + d_dy1[i]= d_dy0[i]; + d_y0[i] = corrct_y; + d_y[i] = corrct_y; + }); + Kokkos::fence(); + + return true; +} \ No newline at end of file diff --git a/src/Integration/AdamsMoulton4/AdamsMoulton4.H b/src/Integration/AdamsMoulton4/AdamsMoulton4.H new file mode 100644 index 00000000..b54750f1 --- /dev/null +++ b/src/Integration/AdamsMoulton4/AdamsMoulton4.H @@ -0,0 +1,185 @@ +/*------------------------------- 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 __AdamsMoulton4_H__ +#define __AdamsMoulton4_H__ + + +#include "integration.H" +#include "pointFields.H" + +namespace pFlow +{ + + +class AdamsMoulton4 +: + public integration +{ +protected: + + realx3PointField_D& y0_; + + realx3PointField_D& dy0_; + + realx3PointField_D& dy1_; + + realx3PointField_D& dy2_; + + using rpIntegration = Kokkos::RangePolicy< + DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType + >; +public: + + // type info + TypeName("AdamsMoulton4"); + + //// - Constructors + AdamsMoulton4( + const word& baseName, + repository& owner, + const pointStructure& pStruct, + const word& method); + + virtual ~AdamsMoulton4()=default; + + // - add a virtual constructor + add_vCtor( + integration, + AdamsMoulton4, + word); + + + //// - Methods + bool predict(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + + bool correct(real dt, realx3Vector_D& y, realx3Vector_D& dy) override; + + bool setInitialVals( + const int32IndexContainer& newIndices, + const realx3Vector& y) override; + + bool needSetInitialVals()const override + { + return true; + } + + uniquePtr clone()const override + { + return makeUnique(*this); + } + + bool predictAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + + template + bool predictRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP); + + bool intAll(real dt, realx3Vector_D& y, realx3Vector_D& dy, range activeRng); + + template + bool intRange(real dt, realx3Vector_D& y, realx3Vector_D& dy, activeFunctor activeP ); + +}; + + +template +bool AdamsMoulton4::predictRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ) +{ + auto d_dy = dy.deviceVectorAll(); + auto d_y = y.deviceVectorAll(); + + auto d_y0 = y0_.deviceVectorAll(); + auto d_dy0 = dy0_.deviceVectorAll(); + auto d_dy1 = dy1_.deviceVectorAll(); + auto d_dy2 = dy2_.deviceVectorAll(); + + auto activeRng = activeP.activeRange(); + + Kokkos::parallel_for( + "AdamsMoulton4::predictRange", + rpIntegration (activeRng.first, activeRng.second), + LAMBDA_HD(int32 i){ + if(activeP(i)) + { + d_dy0[i] = d_dy[i]; + d_y[i] = d_y0[i] + dt*( + static_cast(23.0 /12.0 ) * d_dy[i] + - static_cast(16.0 / 12.0) * d_dy1[i] + + static_cast( 5.0 / 12.0) * d_dy2[i]); + } + }); + Kokkos::fence(); + + return true; + +} + + +template +bool pFlow::AdamsMoulton4::intRange( + real dt, + realx3Vector_D& y, + realx3Vector_D& dy, + activeFunctor activeP ) +{ + + auto d_dy = dy.deviceVectorAll(); + auto d_y = y.deviceVectorAll(); + + auto d_dy0 = dy0_.deviceVectorAll(); + auto d_y0 = y0_.deviceVectorAll(); + auto d_dy1 = dy1_.deviceVectorAll(); + auto d_dy2 = dy2_.deviceVectorAll(); + + auto activeRng = activeP.activeRange(); + + Kokkos::parallel_for( + "AdamsMoulton4::correct", + rpIntegration (activeRng.first, activeRng.second), + LAMBDA_HD(int32 i){ + if( activeP(i)) + { + auto corrct_y = d_y0[i] + dt*( + static_cast(9.0/24.0)*d_dy[i] + + static_cast(19.0/24.0)*d_dy0[i] + - static_cast( 5.0/24.0)*d_dy1[i] + + static_cast( 1.0/24.0)*d_dy2[i]); + + d_dy2[i]= d_dy1[i]; + d_dy1[i]= d_dy0[i]; + d_y0[i] = corrct_y; + d_y[i] = corrct_y; + } + }); + Kokkos::fence(); + + + return true; +} + +} // pFlow + +#endif //__integration_H__ diff --git a/src/Integration/CMakeLists.txt b/src/Integration/CMakeLists.txt index 183b0477..63088123 100644 --- a/src/Integration/CMakeLists.txt +++ b/src/Integration/CMakeLists.txt @@ -6,6 +6,7 @@ AdamsBashforth4/AdamsBashforth4.C AdamsBashforth3/AdamsBashforth3.C AdamsBashforth2/AdamsBashforth2.C AdamsMoulton3/AdamsMoulton3.C +AdamsMoulton4/AdamsMoulton4.C ) set(link_libs Kokkos::kokkos phasicFlow)