diff --git a/src/TNL/Containers/Algorithms/CommonVectorOperations.hpp b/src/TNL/Containers/Algorithms/CommonVectorOperations.hpp index a9e9161cf1f7eb9a2d885a1aaca9a27e8841e549..b8646af9e68d5d4c519549e408c73762eadc6766 100644 --- a/src/TNL/Containers/Algorithms/CommonVectorOperations.hpp +++ b/src/TNL/Containers/Algorithms/CommonVectorOperations.hpp @@ -48,7 +48,7 @@ getVectorMin( const Vector& v ) using IndexType = typename Vector::IndexType; const auto* data = v.getData(); - auto fetch = [=] __cuda_callable__ ( IndexType i ) { return data[ i ]; }; + auto fetch = [=] __cuda_callable__ ( IndexType i ) -> RealType { return data[ i ]; }; auto reduction = [=] __cuda_callable__ ( ResultType& a, const ResultType& b ) { a = TNL::min( a, b ); }; auto volatileReduction = [=] __cuda_callable__ ( volatile ResultType& a, volatile ResultType& b ) { a = TNL::min( a, b ); }; return Reduction< DeviceType >::reduce( v.getSize(), reduction, volatileReduction, fetch, std::numeric_limits< ResultType >::max() ); diff --git a/src/TNL/Solvers/IterativeSolver.h b/src/TNL/Solvers/IterativeSolver.h index 201b20ae6398a5219aa37e6809a10a538f41a84a..be0ee9f1d2558b4481a52c71a21e3d785e729c24 100644 --- a/src/TNL/Solvers/IterativeSolver.h +++ b/src/TNL/Solvers/IterativeSolver.h @@ -15,7 +15,7 @@ #include <TNL/Solvers/IterativeSolverMonitor.h> namespace TNL { -namespace Solvers { +namespace Solvers { template< typename Real, typename Index > class IterativeSolver diff --git a/src/TNL/Solvers/ODE/Euler.h b/src/TNL/Solvers/ODE/Euler.h index 5971fb8c040d261d8a8fde62b08392ed30d7882a..c782b08cd96226d353d0ca30eb9bbe4ff7c8ad81 100644 --- a/src/TNL/Solvers/ODE/Euler.h +++ b/src/TNL/Solvers/ODE/Euler.h @@ -25,13 +25,13 @@ class Euler : public ExplicitSolver< Problem > { public: - typedef Problem ProblemType; - typedef typename Problem :: DofVectorType DofVectorType; - typedef typename Problem :: RealType RealType; - typedef typename Problem :: DeviceType DeviceType; - typedef typename Problem :: IndexType IndexType; - typedef Pointers::SharedPointer< DofVectorType, DeviceType > DofVectorPointer; - + using ProblemType = Problem; + using DofVectorType = typename ProblemType::DofVectorType; + using RealType = typename ProblemType::RealType; + using DeviceType = typename ProblemType::DeviceType; + using IndexType = typename ProblemType::IndexType; + using DofVectorView = typename ViewTypeGetter< DofVectorType >::Type; + using DofVectorPointer = Pointers::SharedPointer< DofVectorType, DeviceType >; Euler(); @@ -54,16 +54,13 @@ class Euler : public ExplicitSolver< Problem > RealType tau, RealType& currentResidue ); - - DofVectorPointer k1; + DofVectorPointer _k1; RealType cflCondition; - - //Timer timer, updateTimer; }; } // namespace ODE } // namespace Solvers } // namespace TNL -#include <TNL/Solvers/ODE/Euler_impl.h> +#include <TNL/Solvers/ODE/Euler.hpp> diff --git a/src/TNL/Solvers/ODE/Euler.hpp b/src/TNL/Solvers/ODE/Euler.hpp index 0b9eed1f8e9e7a356ae78d21469cde49bb88c6b1..5e2c4979e7cdc689f221dea1344113670b801e2b 100644 --- a/src/TNL/Solvers/ODE/Euler.hpp +++ b/src/TNL/Solvers/ODE/Euler.hpp @@ -72,14 +72,16 @@ const typename Problem :: RealType& Euler< Problem > :: getCFLCondition() const } template< typename Problem > -bool Euler< Problem > :: solve( DofVectorPointer& u ) +bool Euler< Problem > :: solve( DofVectorPointer& _u ) { /**** * First setup the supporting meshes k1...k5 and k_tmp. */ //timer.start(); - k1->setLike( *u ); - k1->setValue( 0.0 ); + _k1->setLike( *_u ); + DofVectorView k1 = _k1->getView(); + DofVectorView u = _u->getView(); + k1 = 0.0; /**** @@ -103,14 +105,14 @@ bool Euler< Problem > :: solve( DofVectorPointer& u ) * Compute the RHS */ //timer.stop(); - this->problem->getExplicitUpdate( time, currentTau, u, k1 ); + this->problem->getExplicitUpdate( time, currentTau, _u, _k1 ); //timer.start(); RealType lastResidue = this->getResidue(); RealType maxResidue( 0.0 ); if( this -> cflCondition != 0.0 ) { - maxResidue = k1->absMax(); + maxResidue = max( abs( k1 ) ); //k1->absMax(); if( currentTau * maxResidue > this->cflCondition ) { currentTau *= 0.9; @@ -118,7 +120,15 @@ bool Euler< Problem > :: solve( DofVectorPointer& u ) } } RealType newResidue( 0.0 ); - computeNewTimeLevel( u, currentTau, newResidue ); + //computeNewTimeLevel( u, currentTau, newResidue ); + /* auto fetch = [ currentTau, k1, u ] __cuda_callable__ ( IndexType i ) -> RealType { + const RealType add = currentTau * k1[ i ]; + u[ i ] += add; + return TNL::abs( add ); }; + auto reduction = [=] __cuda_callable__ ( RealType& a , const RealType& b ) { a += b; }; + auto volatileReduction = [=] __cuda_callable__ ( volatile RealType& a , const volatile RealType& b ) { a += b; }; + return Containers::Algorithms::Reduction< DeviceType >::reduce( u.getSize(), reduction, volatileReduction, fetch, 0.0 );*/ + this->setResidue( newResidue ); /**** @@ -127,7 +137,7 @@ bool Euler< Problem > :: solve( DofVectorPointer& u ) */ if( currentTau + time == this -> stopTime ) this->setResidue( lastResidue ); time += currentTau; - this->problem->applyBoundaryConditions( time, u ); + this->problem->applyBoundaryConditions( time, _u ); if( ! this->nextIteration() ) return this->checkConvergence(); @@ -159,7 +169,7 @@ bool Euler< Problem > :: solve( DofVectorPointer& u ) } }; -template< typename Problem > +/*template< typename Problem > void Euler< Problem > :: computeNewTimeLevel( DofVectorPointer& u, RealType tau, RealType& currentResidue ) @@ -262,6 +272,7 @@ __global__ void updateUEuler( const Index size, n ); } #endif +*/ } // namespace ODE } // namespace Solvers