diff --git a/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp b/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp index 09733eb4..54af527c 100644 --- a/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp +++ b/src/Integration/AdamsBashforth2/AdamsBashforth2.cpp @@ -50,7 +50,7 @@ bool intAllActive( "AdamsBashforth2::correct", rpIntegration (activeRng.start(), activeRng.end()), LAMBDA_HD(uint32 i){ - d_y[i] = damping*(d_dy[i] + dt*(static_cast(1.5) * d_dy[i] - static_cast(0.5) * d_dy1[i])); + d_y[i] = damping*(d_y[i] + dt*(static_cast(1.5) * d_dy[i] - static_cast(0.5) * d_dy1[i])); d_dy1[i] = d_dy[i]; }); Kokkos::fence(); diff --git a/src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp b/src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp index 25180208..38566c23 100755 --- a/src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp +++ b/src/Interaction/Models/contactForce/cGAbsoluteLinearCF.hpp @@ -25,12 +25,6 @@ Licence: #include "symArrays.hpp" - - - - - - namespace pFlow::cfModels { @@ -159,7 +153,7 @@ protected: { etha_t[i] = -2.0*log( et[i])*sqrt(kt[i]*2/7) / - sqrt(log(pow(et[i],2.0))+ pow(Pi,2.0)); + sqrt(log(pow(et[i],2))+ pow(Pi,2)); } Vector prop("prop", nElem); @@ -285,8 +279,8 @@ public: 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 mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j]; real sqrt_meff = sqrt((mi*mj)/(mi+mj)); @@ -299,29 +293,24 @@ public: 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)))) ) )); + prop.en_ = exp((pow(f_,static_cast(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)); + sqrt(pow(log(prop.en_),2)+ pow(Pi,2)); //REPORT(0)<<"\n en n is : "<(1.5)) * ethan_ * vrn)*Nij; + FCt = ( -pow(f_,3)*prop.kt_ * history.overlap_t_ - sqrt_meff * pow(f_,static_cast(1.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( length(history.overlap_t_) >zero) { if constexpr (limited) { diff --git a/src/Interaction/Models/contactForce/cGNonLinearCF.hpp b/src/Interaction/Models/contactForce/cGNonLinearCF.hpp index a167f8fa..18ae832f 100644 --- a/src/Interaction/Models/contactForce/cGNonLinearCF.hpp +++ b/src/Interaction/Models/contactForce/cGNonLinearCF.hpp @@ -236,7 +236,7 @@ public: const auto prop = nonlinearProperties_(propId_i,propId_j); - const real f = 2.0/( 1.0/cGFi + 1.0/cGFj ); + const real f = 2/( 1/cGFi + 1/cGFj ); real vrn = dot(Vr, Nij); realx3 Vt = Vr - vrn*Nij; @@ -245,7 +245,7 @@ public: real mi = 3*Pi/4*pow(Ri,static_cast(3))*rho_[propId_i]; real mj = 3*Pi/4*pow(Rj,static_cast(3))*rho_[propId_j]; - real Reff = 1.0/(1/Ri + 1/Rj); + real Reff = 1/(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); @@ -258,13 +258,13 @@ public: 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)))) ) )); + en = exp((pow(f,static_cast(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(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)); + sqrt(pow(log(en),2)+ pow(Pi,2)); FCn = ( - Kn*ovrlp_n - sqrt_meff_K_hertz*ethan_*pow(ovrlp_n,static_cast(0.25))*vrn)*Nij; diff --git a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp index d542e70c..af51755d 100755 --- a/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp +++ b/src/Interaction/Models/contactForce/cGRelativeLinearCF.hpp @@ -255,8 +255,8 @@ public: 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 mi = 3*Pi/4*pow(Ri,3)*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,3)*rho_[propId_j]; real sqrt_meff = sqrt((mi*mj)/(mi+mj)); @@ -268,14 +268,26 @@ public: 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)))) ) )); + en = exp( + ( + pow(f_,static_cast(1.5))* + log(prop.en_)* + sqrt( + (1-((pow(log(prop.en_),static_cast(2.0)) + )/ + ( + pow(log(prop.en_),static_cast(2.0))+ + pow(Pi,static_cast(2.0)))))/ + (1-(pow(f_,3)*(pow(log(prop.en_),2))/ + (pow(log(prop.en_),static_cast(2.0))+ + pow(Pi,static_cast(2.0))))) ) )); } real ethan_ = -2.0*log(en)*sqrt(prop.kn_)/ - sqrt(pow(log(en),2.0)+ pow(Pi,2.0)); + sqrt(pow(log(en),2)+ pow(Pi,2)); - FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,0.5) * ethan_ * vrn)*Nij; + FCn = ( -f_*prop.kn_ * ovrlp_n - sqrt_meff * pow(f_,half) * ethan_ * vrn)*Nij; FCt = ( -f_*prop.kt_ * history.overlap_t_); diff --git a/src/Interaction/Models/contactForce/linearCF.hpp b/src/Interaction/Models/contactForce/linearCF.hpp index 88bf2108..bc3a030e 100644 --- a/src/Interaction/Models/contactForce/linearCF.hpp +++ b/src/Interaction/Models/contactForce/linearCF.hpp @@ -139,10 +139,10 @@ protected: ForAll(i , kn) { etha_n[i] = -2.0*log(en[i])*sqrt(kn[i])/ - sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0)); + sqrt(pow(log(en[i]),static_cast(2.0))+ pow(Pi,static_cast(2.0))); etha_t[i] = -2.0*log( et[i]*sqrt(kt[i]) )/ - sqrt(pow(log(et[i]),2.0)+ pow(Pi,2.0)); + sqrt(pow(log(et[i]),static_cast(2.0))+ pow(Pi,static_cast(2.0))); } Vector prop("prop", nElem); @@ -243,8 +243,8 @@ public: 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 mi = 3*Pi/4*pow(Ri,static_cast(3.0))*rho_[propId_i]; + real mj = 3*Pi/4*pow(Rj,static_cast(3.0))*rho_[propId_j]; real sqrt_meff = sqrt((mi*mj)/(mi+mj)); diff --git a/src/Interaction/Models/contactForce/nonLinearCF.hpp b/src/Interaction/Models/contactForce/nonLinearCF.hpp index 3ba1c5c5..4fdd4498 100644 --- a/src/Interaction/Models/contactForce/nonLinearCF.hpp +++ b/src/Interaction/Models/contactForce/nonLinearCF.hpp @@ -131,7 +131,7 @@ protected: // we take out sqrt(meff*K_hertz) here and then consider this term // when calculating damping part. etha_n[i] = -2.2664*log(en[i])/ - sqrt(pow(log(en[i]),2.0)+ pow(Pi,2.0)); + sqrt(pow(log(en[i]),static_cast(2.0))+ pow(Pi,static_cast(2.0))); // no damping for tangential part @@ -255,7 +255,7 @@ public: // apply friction if(ft > ft_fric) { - if( length(history.overlap_t_) >0.0) + if( length(history.overlap_t_) >zero) { if constexpr (limited) { diff --git a/src/Interaction/Models/rolling/grainRolling.hpp b/src/Interaction/Models/rolling/grainRolling.hpp index b8eff130..2c5c31a2 100755 --- a/src/Interaction/Models/rolling/grainRolling.hpp +++ b/src/Interaction/Models/rolling/grainRolling.hpp @@ -96,10 +96,10 @@ public: realx3 w_hat = wi-wj; real w_hat_mag = length(w_hat); - if( !equal(w_hat_mag,0.0) ) + if( !equal(w_hat_mag,zero) ) w_hat /= w_hat_mag; else - w_hat = 0.0; + w_hat = zero; auto Reff = (Ri*Rj)/(Ri+Rj); diff --git a/src/Interaction/Models/rolling/normalRolling.hpp b/src/Interaction/Models/rolling/normalRolling.hpp index 80b8e321..7f9fc026 100644 --- a/src/Interaction/Models/rolling/normalRolling.hpp +++ b/src/Interaction/Models/rolling/normalRolling.hpp @@ -94,10 +94,10 @@ public: realx3 w_hat = wi-wj; real w_hat_mag = length(w_hat); - if( !equal(w_hat_mag,0.0) ) + if( !equal(w_hat_mag,zero) ) w_hat /= w_hat_mag; else - w_hat = 0.0; + w_hat = zero; auto Reff = (Ri*Rj)/(Ri+Rj); diff --git a/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp b/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp index acfad226..e6aabb9b 100644 --- a/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp +++ b/src/Interaction/contactSearch/boundaries/periodicBoundaryContactSearch/wallBoundaryContactSearch.cpp @@ -11,7 +11,7 @@ pFlow::wallBoundaryContactSearch::wallBoundaryContactSearch const ViewType1D &normals ) : - cellExtent_( max(cellExtent, 0.5 ) ), + cellExtent_( max(cellExtent, half ) ), numElements_(numElements), numPoints_(numPoints), vertices_(vertices), diff --git a/src/Interaction/contactSearch/contactSearch/contactSearch.cpp b/src/Interaction/contactSearch/contactSearch/contactSearch.cpp index 7394e295..b384b0bc 100644 --- a/src/Interaction/contactSearch/contactSearch/contactSearch.cpp +++ b/src/Interaction/contactSearch/contactSearch/contactSearch.cpp @@ -31,6 +31,7 @@ pFlow::contactSearch::contactSearch( Timers& timers) : extendedDomainBox_(extDomain), + updateInterval_(dict.getValMax("updateInterval", 1u)), particles_(prtcl), geometry_(geom), bTimer_("Boundary particles contact search", &timers), diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp index 733313b4..916cb117 100644 --- a/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/NBS.cpp @@ -44,8 +44,8 @@ pFlow::NBS::NBS position, flags, diam), - sizeRatio_(max(dict.getVal("sizeRatio"), 1.0)), - cellExtent_(max(dict.getVal("cellExtent"), 0.5)), + sizeRatio_(max(dict.getVal("sizeRatio"), one)), + cellExtent_(max(dict.getVal("cellExtent"), half)), adjustableBox_(dict.getVal("adjustableBox")), NBSLevel0_ ( diff --git a/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp index e1b6e3a7..c565bf87 100644 --- a/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp +++ b/src/Interaction/contactSearch/methods/cellBased/NBS/cellsWallLevel0.cpp @@ -31,7 +31,7 @@ pFlow::cellsWallLevel0::cellsWallLevel0 const ViewType1D& normals ) : - cellExtent_( max(cellExtent, 0.5 ) ), + cellExtent_( max(cellExtent, half ) ), numElements_(numElements), numPoints_(numPoints), vertices_(vertices), diff --git a/src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp b/src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp index d5aa3ca2..1b2a4c58 100644 --- a/src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp +++ b/src/MotionModel/entities/rotatingAxis/rotatingAxis.cpp @@ -47,7 +47,7 @@ pFlow::rotatingAxis::rotatingAxis timeInterval(), line(p1, p2), omega_(omega), - rotating_(!equal(omega,0.0)) + rotating_(!equal(omega,zero)) { } @@ -58,7 +58,7 @@ pFlow::real pFlow::rotatingAxis::setOmega(real omega) { auto tmp = omega_; omega_ = omega; - rotating_ = !equal(omega, 0.0); + rotating_ = !equal(omega, zero); return tmp; } @@ -72,7 +72,7 @@ bool pFlow::rotatingAxis::read if(!timeInterval::read(dict))return false; if(!line::read(dict)) return false; - real omega = dict.getValOrSet("omega", 0.0); + real omega = dict.getValOrSet("omega", zero); setOmega(omega); return true; diff --git a/src/phasicFlow/containers/symArrayHD/symArrayHD.hpp b/src/phasicFlow/containers/symArrayHD/symArrayHD.hpp index 9c7d1a67..73cc4dbf 100644 --- a/src/phasicFlow/containers/symArrayHD/symArrayHD.hpp +++ b/src/phasicFlow/containers/symArrayHD/symArrayHD.hpp @@ -239,7 +239,7 @@ public: { real nr = 0.5*(sqrt(8.0*nElem+1)-1); n = static_cast(nr); - if( equal(nr-static_cast(n) , 0.0) ) return true; + if( equal(nr-static_cast(n) , zero) ) return true; return false; } diff --git a/src/phasicFlow/ranges/stridedRange.hpp b/src/phasicFlow/ranges/stridedRange.hpp index 02df31d2..89d6f684 100644 --- a/src/phasicFlow/ranges/stridedRange.hpp +++ b/src/phasicFlow/ranges/stridedRange.hpp @@ -111,7 +111,7 @@ public: { if(!isInRange(val)) return false; - if(T dist = val-begin_; abs(dist%stride_)<= epsilon) return true; + if(const T dist = val-begin_; abs(dist%stride_)<=epsilon) return true; if(equal(val,begin_))return true; if(equal(val,end_))return true; return false; @@ -151,19 +151,34 @@ public: template<> inline -bool stridedRange::isMember(real val, real epsilon)const +bool stridedRange::isMember(float val, float epsilon)const { - if(!isInRange(val)) return false; + /*if(!isInRange(val)) return false; real dist = val-begin_; if(abs( (dist-(static_cast((dist+0.01*epsilon)/stride_)*stride_)) )<= epsilon) return true; if(equal(val,begin_))return true; if(equal(val,end_))return true; + return false;*/ + if(!isInRange(val)) return false; + if(const float dist = val-begin_; abs(remainder(dist,stride_)<= epsilon)) return true; + if(equal(val,begin_))return true; + if(equal(val,end_))return true; return false; } +template<> +inline +bool stridedRange::isMember(double val, double epsilon)const +{ + if(!isInRange(val)) return false; + if(const double dist = val-begin_; abs(remainder(dist,stride_)<= epsilon)) return true; + if(equal(val,begin_))return true; + if(equal(val,end_))return true; + return false; +} } diff --git a/src/phasicFlow/repository/Time/timeControl.cpp b/src/phasicFlow/repository/Time/timeControl.cpp index 98cd24ce..2990b768 100644 --- a/src/phasicFlow/repository/Time/timeControl.cpp +++ b/src/phasicFlow/repository/Time/timeControl.cpp @@ -25,7 +25,7 @@ Licence: bool pFlow::timeControl::screenReport()const { - return screenReportInterval_.isMember(currentIter_); + return screenReportInterval_.isMember(ti_.currentIter()); } pFlow::timeControl::timeControl @@ -36,16 +36,16 @@ pFlow::timeControl::timeControl ti_(dict), startTime_ ( - dict.getVal("startTime") + dict.getVal("startTime") ), endTime_ ( - dict.getVal("endTime") + dict.getVal("endTime") ), stopAt_(endTime_), saveInterval_ ( - dict.getVal("saveInterval") + dict.getVal("saveInterval") ), lastSaved_(startTime_), performSorting_ @@ -65,50 +65,40 @@ pFlow::timeControl::timeControl pFlow::timeControl::timeControl( dictionary& dict, - real startTime, - real endTime, - real saveInterval, + timeValue startTime, + timeValue endTime, + timeValue saveInterval, word startTimeName) : - dt_ - ( - dict.getVal("dt") - ), + ti_(startTime, dict), startTime_(startTime), endTime_(endTime), stopAt_(endTime_), - currentTime_(startTime_), saveInterval_(saveInterval), - lastSaved_(startTime_), - currentIter_(0), - timePrecision_ - ( - dict.getValOrSet("timePrecision", 4) - ), + lastSaved_(startTime_), managedExternaly_(true), - timeName_(startTimeName), - timersReportInterval_ - ( - startTime_, - dict.getValOrSet("timersReportInterval", 0.04) - ), performSorting_ ( dict.getValOrSet("performSorting", Logical("No")) - ), - sortingInterval_ - ( - startTime_, - dict.getValOrSet("sortingInterval", static_cast(1.0)) ) { + timeName_ = startTimeName; + + timersReportInterval_ = timeStridedRange( + startTime_, + dict.getValOrSet("timersReportInterval", 0.04)); + + sortingInterval_ = timeStridedRange( + startTime_, + dict.getValOrSet("sortingInterval", 1.0)); + checkForOutputToFile(); } -pFlow::real pFlow::timeControl::setTime(real t) +pFlow::timeValue pFlow::timeControl::setTime(timeValue t) { - real tmp = currentTime_; - currentTime_ = t; + timeValue tmp = ti_.currentTime(); + ti_.currentTime_ = t; lastSaved_ = t; checkForOutputToFile(); return tmp; @@ -124,15 +114,15 @@ pFlow::word pFlow::timeControl::timeName()const bool pFlow::timeControl::finalTime()const { - if( currentTime_ >= endTime_ ) return true; - if( std::abs(currentTime_-endTime_) < 0.5*dt_ )return true; + if( ti_.currentTime_ >= endTime_ ) return true; + if( std::abs(ti_.currentTime_-endTime_) < 0.5*ti_.dt_ )return true; return false; } bool pFlow::timeControl::reachedStopAt()const { - if( currentTime_ >= stopAt_ ) return true; - if( std::abs(currentTime_-stopAt_) < 0.5*dt_ )return true; + if( ti_.currentTime_ >= stopAt_ ) return true; + if( std::abs(ti_.currentTime_-stopAt_) < 0.5*ti_.dt_ )return true; return false; } @@ -142,22 +132,22 @@ void pFlow::timeControl::checkForOutputToFile() bool save = false; if(managedExternaly_) { - if( std::abs(currentTime_-writeTime_) < 0.5*dt_) + if( std::abs(ti_.currentTime_-writeTime_) < 0.5*ti_.dt_) { save = true; - lastSaved_ = currentTime_; + lastSaved_ = ti_.currentTime_; } } else { - if ( std::abs(currentTime_ - lastSaved_ - saveInterval_) < 0.5 * dt_ ) + if ( std::abs(ti_.currentTime_ - lastSaved_ - saveInterval_) < 0.5 * ti_.dt_ ) { - lastSaved_ = currentTime_; + lastSaved_ = ti_.currentTime_; save = true; } - else if( std::abs(currentTime_ - lastSaved_) < std::min( pow(10.0,-1.0*timePrecision_), static_cast(0.5 *dt_)) ) + else if( std::abs(ti_.currentTime_ - lastSaved_) < std::min( std::pow(10.0,-1.0*ti_.precision()), 0.5 *ti_.dt_) ) { - lastSaved_ = currentTime_; + lastSaved_ = ti_.currentTime_; save = true; } @@ -168,13 +158,13 @@ void pFlow::timeControl::checkForOutputToFile() bool pFlow::timeControl::timersReportTime()const { - if(currentIter_<=1)return false; - return timersReportInterval_.isMember(currentTime_, 0.55*dt_); + if(ti_.currentIter_<=1)return false; + return timersReportInterval_.isMember(ti_.currentTime_, 0.55*ti_.dt_); } bool pFlow::timeControl::sortTime()const { - return performSorting_()&&sortingInterval_.isMember(currentTime_,dt_); + return performSorting_()&&sortingInterval_.isMember(ti_.currentTime_,ti_.dt_); } void pFlow::timeControl::setSaveTimeFolder( @@ -192,10 +182,10 @@ bool pFlow::timeControl::operator ++(int) { if( reachedStopAt() ) return false; + // increament iteration number - currentIter_++; + ti_.martchDT(); - currentTime_ += dt_; if(screenReport() && !managedExternaly_) { REPORT(0)<<"Time (s): "<; private: //// - Data members - // - integration time step - real dt_; - - // - start time of simulation - real startTime_; + // - time information + timeInfo ti_; + + // - start time of simulation + timeValue startTime_; - // - end time of simulation - real endTime_; + // - end time of simulation + timeValue endTime_; - // - stopAt - real stopAt_; + // - stopAt + timeValue stopAt_; - // - current time of simulation - real currentTime_; + // - time interval for time folder output + timeValue saveInterval_; - // - time interval for time folder output - real saveInterval_; + // - the last time folder that was saved + timeValue lastSaved_; - // - the last time folder that was saved - real lastSaved_; + bool managedExternaly_ = false; - // - current iteration number (for this execution) - int32 currentIter_; + bool outputToFile_ = false; - // - time precision for output time folders - int32 timePrecision_; + Logical performSorting_; - bool managedExternaly_ = false; + static + inline int32StridedRagne screenReportInterval_ ={0,100}; - word timeName_ = "wrongSettings"; // for managedExternamly + static + inline timeStridedRange timersReportInterval_{0,100}; - real writeTime_ = 0; // for managedExternamly + static + inline timeStridedRange sortingInterval_{0,100}; + + static + inline word timeName_ = "wrongSettings"; // for managedExternamly - realStridedRange timersReportInterval_; + static + inline timeValue writeTime_ = 0; // for managedExternamly - Logical performSorting_; - - realStridedRange sortingInterval_; - - int32StridedRagne screenReportInterval_ ={0,100}; - - bool outputToFile_ = false; - - void checkForOutputToFile(); - - bool screenReport()const; + void checkForOutputToFile(); + + bool screenReport()const; public: @@ -94,22 +92,22 @@ public: timeControl( dictionary& dict, - real startTime, - real endTime, - real saveInterval, + timeValue startTime, + timeValue endTime, + timeValue saveInterval, word startTimeName); virtual ~timeControl() = default; - real dt()const + timeValue dt()const { - return dt_; + return ti_.dt(); } - real setTime(real t); + timeValue setTime(timeValue t); - void setStopAt(real sT) + void setStopAt(timeValue sT) { if(managedExternaly_) { @@ -117,37 +115,27 @@ public: } } - real startTime()const + timeValue startTime()const { return startTime_; } word timeName()const; - real currentTime() const + timeValue currentTime() const { - return currentTime_; + return ti_.currentTime(); } word currentTimeWord(bool forSave = true)const { - return real2FixedStripZeros( currentTime(), timePrecision()); - /*if(forSave) - { - if(!managedExternaly_) - - else - return timeName_; - } - else - { - return real2FixedStripZeros( currentTime(), timePrecision()); - }*/ + return ti_.timeName(); + } int32 currentIter()const { - return currentIter_; + return ti_.currentIter(); } bool finalTime()const; @@ -179,17 +167,17 @@ public: bool saveToFile, const word& timeName = "wrongTimeFolder"); - int32 timePrecision()const - { - return timePrecision_; - } - inline - timeInfo TimeInfo()const + const timeInfo& TimeInfo()const { - return {static_cast(currentIter_), currentTime_, dt_}; + return ti_; } + static + int32 timePrecision() + { + return timeInfo::precision(); + } }; diff --git a/src/phasicFlow/repository/Time/timeInfo.hpp b/src/phasicFlow/repository/Time/timeInfo.hpp index 34e30d90..c983d577 100644 --- a/src/phasicFlow/repository/Time/timeInfo.hpp +++ b/src/phasicFlow/repository/Time/timeInfo.hpp @@ -21,6 +21,7 @@ Licence: #define __timeInfo__hpp_ #include "types.hpp" +#include "dictionary.hpp" namespace pFlow { @@ -28,34 +29,88 @@ namespace pFlow class timeInfo { private: + + friend class timeControl; - uint32 currentIter_; + // - current iteration number (for this execution) + int32 currentIter_; - real currentTime_; + // - current time of simulation + timeValue currentTime_; - real dt_; + // - integration time step + timeValue dt_; + + inline static int32 presicion_ = 5; public: - timeInfo(uint32 cIter, real cTime, real dt) + timeInfo(int32 cIter, timeValue cTime, timeValue dt) : currentIter_(cIter), currentTime_(cTime), dt_(dt) { } - inline const real& t() const + timeInfo(const dictionary& dict) + : + currentIter_(0), + currentTime_(dict.getVal("startTime")), + dt_( dict.getVal("dt")) + { + presicion_ = dict.getValOrSet("timePrecision",5); + } + + timeInfo(timeValue currentTime, const dictionary& dict) + : + currentIter_(0), + currentTime_(currentTime), + dt_( dict.getVal("dt")) + { + presicion_ = dict.getValOrSet("timePrecision",5); + } + + inline const timeValue& currentTime()const { return currentTime_; } - inline const real& dt() const + + inline const timeValue& t() const + { + return currentTime_; + } + inline const timeValue& dt() const { return dt_; } - inline const uint32& iter() const + inline const int32& iter() const { return currentIter_; } + + inline const int32& currentIter() const + { + return currentIter_; + } + + inline + void martchDT() + { + currentIter_++; + currentTime_ += dt_; + } + + inline + word timeName()const + { + return real2FixedStripZeros(currentTime_, presicion_); + } + + static + int32 precision() + { + return presicion_; + } }; } // namespace pFlow diff --git a/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp b/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp index 89e497fb..cc55eab7 100644 --- a/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp +++ b/src/phasicFlow/structuredData/infinitePlane/infinitePlane.cpp @@ -29,7 +29,7 @@ pFlow::infinitePlane::infinitePlane { auto ln = cross(p2-p1, p3-p1); - if( equal(ln.length(),0.0) ) + if( equal(ln.length(),zero) ) { fatalErrorInFunction<< "invalid input to form a infinte wall "<< realx3x3(p1,p2,p3)<= 0.0 && d <= dist; + return d >= zero && d <= dist; } INLINE_FUNCTION_HD bool inNegativeDistance(const realx3& p, real dist)const { real d = pointFromPlane(p); - return d < 0.0 && d >= -dist; + return d < zero && d >= -dist; } INLINE_FUNCTION_HD bool pointOnPlane(const realx3& p)const { - return equal(pointFromPlane(p),0.0); + return equal(pointFromPlane(p),zero); } INLINE_FUNCTION_HD diff --git a/src/phasicFlow/triSurface/triangleFunctions.hpp b/src/phasicFlow/triSurface/triangleFunctions.hpp index c6812d2f..2a4cc87f 100644 --- a/src/phasicFlow/triSurface/triangleFunctions.hpp +++ b/src/phasicFlow/triSurface/triangleFunctions.hpp @@ -38,8 +38,8 @@ INLINE_FUNCTION_HD realx3 normal(const realx3& p1, const realx3& p2, const realx3& p3) { auto n = cross(p2-p1, p3-p1); - if( equal(n.length(), 0.0) ) - return realx3(0); + if( equal(n.length(), zero) ) + return zero3; else return normalize(n); } @@ -52,11 +52,9 @@ bool valid const realx3& p3 ) { - return !equal(cross(p2-p1, p3-p1).length(), 0.0); + return !equal(cross(p2-p1, p3-p1).length(), zero); } - - } #endif diff --git a/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp b/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp index dedf002d..09835aa7 100755 --- a/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp +++ b/src/phasicFlow/types/basicTypes/bTypesFunctions.cpp @@ -81,7 +81,7 @@ pFlow::isNo(const word& str) } pFlow::word -pFlow::real2Fixed(const real& v, int32 numPrecision) +pFlow::real2Fixed(const float& v, int32 numPrecision) { std::stringstream ss; @@ -90,7 +90,32 @@ pFlow::real2Fixed(const real& v, int32 numPrecision) } pFlow::word -pFlow::real2Word(const real& v, int32 numPrecision) +pFlow::real2Fixed(const double& v, int32 numPrecision) +{ + std::stringstream ss; + + ss << std::fixed << std::setprecision(numPrecision) << v; + return ss.str(); +} + +pFlow::word +pFlow::real2Word(const float& v, int32 numPrecision) +{ + std::stringstream ss; + if (abs(v) < verySmallValue) + { + ss << "0"; + } + else + { + ss << std::setprecision(numPrecision) << v; + } + + return ss.str(); +} + +pFlow::word +pFlow::real2Word(const double& v, int32 numPrecision) { std::stringstream ss; if (abs(v) < verySmallValue) @@ -145,7 +170,14 @@ pFlow::removeDecimalZeros(const word& str) } pFlow::word -pFlow::real2FixedStripZeros(const real& v, int32 numPrecision) +pFlow::real2FixedStripZeros(const float& v, int32 numPrecision) +{ + word strVal = real2Fixed(v, numPrecision); + return removeDecimalZeros(strVal); +} + +pFlow::word +pFlow::real2FixedStripZeros(const double& v, int32 numPrecision) { word strVal = real2Fixed(v, numPrecision); return removeDecimalZeros(strVal); diff --git a/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp b/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp index 5ad36fbb..764866d2 100755 --- a/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp +++ b/src/phasicFlow/types/basicTypes/bTypesFunctions.hpp @@ -39,6 +39,8 @@ inline const real zero = 0.0; /// one real variable inline const real one = 1.0; +inline const real half = 0.5; + /// zero int32 variable inline const int32 zero32 = 0; @@ -76,11 +78,17 @@ isNo(const word& str); /// Convert floating point variable to string with fixed number of precisions word -real2Fixed(const real& v, int32 numPrecision = 6); +real2Fixed(const float& v, int32 numPrecision = 6); + +word +real2Fixed(const double& v, int32 numPrecision = 6); /// Convert floating point variable to string with general format word -real2Word(const real& v, int32 numPrecision = 6); +real2Word(const float& v, int32 numPrecision = 6); + +word +real2Word(const double& v, int32 numPrecision = 6); /// Convert int32 to word word @@ -92,7 +100,10 @@ removeDecimalZeros(const word& str); /// Convert to fixed point variable and remove zeros word -real2FixedStripZeros(const real& v, int32 numPrecision = 6); +real2FixedStripZeros(const float& v, int32 numPrecision = 6); + +word +real2FixedStripZeros(const double& v, int32 numPrecision = 6); /// Output word @@ -225,7 +236,14 @@ readValue(const word& w, bool& val) INLINE_FUNCTION_HD bool -equal(const real& s1, const real& s2, real tol = smallValue) +equal(const float& s1, const float& s2, float tol = smallValue) +{ + return abs(s1 - s2) <= tol; +} + +INLINE_FUNCTION_HD +bool +equal(const double& s1, const double& s2, double tol = smallValue) { return abs(s1 - s2) <= tol; } diff --git a/src/phasicFlow/types/basicTypes/builtinTypes.hpp b/src/phasicFlow/types/basicTypes/builtinTypes.hpp index b10827c5..ab8673ed 100755 --- a/src/phasicFlow/types/basicTypes/builtinTypes.hpp +++ b/src/phasicFlow/types/basicTypes/builtinTypes.hpp @@ -63,6 +63,8 @@ using size_t = std::size_t; using word = std::string; +using timeValue = double; + inline const int numBytesForReal__ = sizeof(real); inline word diff --git a/src/phasicFlow/types/basicTypes/math.hpp b/src/phasicFlow/types/basicTypes/math.hpp index fbe3a6f1..1dce56e5 100755 --- a/src/phasicFlow/types/basicTypes/math.hpp +++ b/src/phasicFlow/types/basicTypes/math.hpp @@ -42,8 +42,19 @@ namespace pFlow { INLINE_FUNCTION_HD -real -abs(real x) +double +abs(double x) +{ +#ifdef __CUDACC__ + return ::fabs(x); +#else + return std::fabs(x); +#endif +} + +INLINE_FUNCTION_HD +float +abs(float x) { #ifdef __CUDACC__ return ::fabs(x); @@ -73,8 +84,18 @@ abs(int32 x) #endif } -INLINE_FUNCTION_HD real -mod(real x, real y) +INLINE_FUNCTION_HD float +mod(float x, float y) +{ +#ifdef __CUDACC__ + return ::fmod(x, y); +#else + return std::fmod(x, y); +#endif +} + +INLINE_FUNCTION_HD double +mod(double x, double y) { #ifdef __CUDACC__ return ::fmod(x, y); @@ -108,8 +129,18 @@ mod(uint32 x, uint32 y) return x % y; } -INLINE_FUNCTION_HD real -remainder(real x, real y) +INLINE_FUNCTION_HD float +remainder(float x, float y) +{ +#ifdef __CUDACC__ + return ::remainder(x,y); +#else + return std::remainder(x, y); +#endif +} + +INLINE_FUNCTION_HD double +remainder(double x, double y) { #ifdef __CUDACC__ return ::remainder(x,y); @@ -120,8 +151,19 @@ remainder(real x, real y) INLINE_FUNCTION_HD -real -exp(real x) +float +exp(float x) +{ +#ifdef __CUDACC__ + return ::exp(x); +#else + return std::exp(x); +#endif +} + +INLINE_FUNCTION_HD +double +exp(double x) { #ifdef __CUDACC__ return ::exp(x); @@ -132,8 +174,17 @@ exp(real x) INLINE_FUNCTION_HD -real -log(real x) +float log(float x) +{ +#ifdef __CUDACC__ + return ::log(x); +#else + return std::log(x); +#endif +} + +INLINE_FUNCTION_HD +double log(double x) { #ifdef __CUDACC__ return ::log(x); @@ -144,8 +195,7 @@ log(real x) INLINE_FUNCTION_HD -real -log10(real x) +float log10(float x) { #ifdef __CUDACC__ return ::log10(x); @@ -155,8 +205,27 @@ log10(real x) } INLINE_FUNCTION_HD -real -pow(real x, real y) +double log10(double x) +{ +#ifdef __CUDACC__ + return ::log10(x); +#else + return std::log10(x); +#endif +} + +INLINE_FUNCTION_HD +float pow(float x, float y) +{ +#ifdef __CUDACC__ + return ::pow(x,y); +#else + return std::pow(x, y); +#endif +} + +INLINE_FUNCTION_HD +double pow(double x, double y) { #ifdef __CUDACC__ return ::pow(x,y); @@ -167,8 +236,17 @@ pow(real x, real y) INLINE_FUNCTION_HD -real -sqrt(real x) +float sqrt(float x) +{ +#ifdef __CUDACC__ + return ::sqrt(x); +#else + return std::sqrt(x); +#endif +} + +INLINE_FUNCTION_HD +double sqrt(double x) { #ifdef __CUDACC__ return ::sqrt(x); @@ -180,8 +258,7 @@ sqrt(real x) INLINE_FUNCTION_HD -real -cbrt(real x) +float cbrt(float x) { #ifdef __CUDACC__ return ::cbrt(x); @@ -191,8 +268,17 @@ cbrt(real x) } INLINE_FUNCTION_HD -real -sin(real x) +double cbrt(double x) +{ +#ifdef __CUDACC__ + return ::cbrt(x); +#else + return std::cbrt(x); +#endif +} + +INLINE_FUNCTION_HD +float sin(float x) { #ifdef __CUDACC__ return ::sin(x); @@ -202,8 +288,17 @@ sin(real x) } INLINE_FUNCTION_HD -real -cos(real x) +double sin(double x) +{ +#ifdef __CUDACC__ + return ::sin(x); +#else + return std::sin(x); +#endif +} + +INLINE_FUNCTION_HD +real cos(real x) { #ifdef __CUDACC__ return ::cos(x); @@ -215,8 +310,7 @@ cos(real x) INLINE_FUNCTION_HD -real -tan(real x) +real tan(real x) { #ifdef __CUDACC__ return ::tan(x); @@ -226,8 +320,7 @@ tan(real x) } INLINE_FUNCTION_HD -real -asin(real x) +real asin(real x) { #ifdef __CUDACC__ return ::asin(x); @@ -237,8 +330,7 @@ asin(real x) } INLINE_FUNCTION_HD -real -acos(real x) +real acos(real x) { #ifdef __CUDACC__ return ::acos(x); @@ -248,8 +340,7 @@ acos(real x) } INLINE_FUNCTION_HD -real -atan(real x) +real atan(real x) { #ifdef __CUDACC__ return ::atan(x); @@ -269,8 +360,7 @@ atan2(real y, real x) } INLINE_FUNCTION_HD -real -sinh(real x) +real sinh(real x) { #ifdef __CUDACC__ return ::sinh(x); @@ -280,8 +370,7 @@ sinh(real x) } INLINE_FUNCTION_HD -real -cosh(real x) +real cosh(real x) { #ifdef __CUDACC__ return ::cosh(x); @@ -322,8 +411,7 @@ acosh(real x) } INLINE_FUNCTION_HD -real -atanh(real x) +real atanh(real x) { #ifdef __CUDACC__ return ::atanh(x); @@ -333,8 +421,18 @@ atanh(real x) } -INLINE_FUNCTION_HD real -min(real x, real y) +INLINE_FUNCTION_HD +float min(float x, float y) +{ +#ifdef __CUDACC__ + return ::fmin(x, y); +#else + return std::fmin(x, y); +#endif +} + +INLINE_FUNCTION_HD +double min(double x, double y) { #ifdef __CUDACC__ return ::fmin(x, y); @@ -344,8 +442,7 @@ min(real x, real y) } INLINE_FUNCTION_HD -int64 -min(int32 x, int32 y) +int64 min(int32 x, int32 y) { #ifdef __CUDACC__ return ::min(x,y); @@ -356,8 +453,7 @@ min(int32 x, int32 y) INLINE_FUNCTION_HD -int64 -min(int64 x, int64 y) +int64 min(int64 x, int64 y) { #ifdef __CUDACC__ return ::min(x,y); @@ -369,8 +465,17 @@ min(int64 x, int64 y) INLINE_FUNCTION_HD -uint64 -min(uint64 x, uint64 y) +uint64 min(uint64 x, uint64 y) +{ +#ifdef __CUDACC__ + return ::min(x,y); +#else + return std::min(x, y); +#endif +} + +INLINE_FUNCTION_HD +uint32 min(uint32 x, uint32 y) { #ifdef __CUDACC__ return ::min(x,y); @@ -380,20 +485,18 @@ min(uint64 x, uint64 y) } - -INLINE_FUNCTION_HD uint32 -min(uint32 x, uint32 y) +INLINE_FUNCTION_HD +float max(float x, float y) { #ifdef __CUDACC__ - return ::min(x,y); + return ::fmax(x, y); #else - return std::min(x, y); + return std::fmax(x, y); #endif } - -INLINE_FUNCTION_HD real -max(real x, real y) +INLINE_FUNCTION_HD +double max(double x, double y) { #ifdef __CUDACC__ return ::fmax(x, y); diff --git a/utilities/postprocessPhasicFlow/postprocessPhasicFlow.cpp b/utilities/postprocessPhasicFlow/postprocessPhasicFlow.cpp index c53fafc5..49762836 100644 --- a/utilities/postprocessPhasicFlow/postprocessPhasicFlow.cpp +++ b/utilities/postprocessPhasicFlow/postprocessPhasicFlow.cpp @@ -103,7 +103,7 @@ int main(int argc, char** argv ) continue; } - if( !withZeroFolder && pFlow::equal(folders.time() , 0.0))continue; + if( !withZeroFolder && pFlow::equal(folders.time() , pFlow::zero))continue; post.processTimeFolder(folders);