Commit 7bb01975 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Fixing the semi-implicit solver for CUDA/

parent e08b9c6f

src/matrices/tnlEllpackMatrix_impl.h

+4 −3
Original line number Diff line number Diff line
@@ -636,15 +636,16 @@ void tnlEllpackMatrix< Real, Device, Index >::print( ostream& str ) const 
   for( IndexType row = 0; row < this->getRows(); row++ )

   {

      str <<"Row: " << row << " -> ";

      IndexType i( row * this->rowLengths );

      const IndexType rowEnd( i + this->rowLengths );

      IndexType i = DeviceDependentCode::getRowBegin( *this, row );

      const IndexType rowEnd = DeviceDependentCode::getRowEnd( *this, row );

      const IndexType step = DeviceDependentCode::getElementStep( *this );

      while( i < rowEnd &&

             this->columnIndexes.getElement( i ) < this->columns &&

             this->columnIndexes.getElement( i ) != this->getPaddingIndex() )

      {

         const Index column = this->columnIndexes.getElement( i );

         str << " Col:" << column << "->" << this->values.getElement( i ) << "\t";

         i++;

         i += step;

      }

      str << endl;

   }

src/matrices/tnlMatrixSetter.h

+2 −0
Original line number Diff line number Diff line
@@ -151,6 +151,7 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >, 
                                  const IndexType index,

                                  const CoordinatesType& coordinates )

         {

             //printf("Matrix setter: Index = %d \n", index );

            ( *userData.rowLengths )[ index ] =

                     userData.boundaryConditions->getLinearSystemRowLength( mesh, index, coordinates );

         }
@@ -182,6 +183,7 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >, 
                                  const IndexType index,

                                  const CoordinatesType& coordinates )

         {

            // printf("Matrix setter: Index = %d \n", index );

            ( *userData.rowLengths )[ index ] =

                     userData.differentialOperator->getLinearSystemRowLength( mesh, index, coordinates );

         }

src/mesh/tnlGrid2D.h

+3 −0
Original line number Diff line number Diff line
@@ -302,6 +302,9 @@ class tnlGrid< 2, Real, Device, Index > : public tnlObject 


   protected:



#ifdef HAVE_CUDA

   __device__ __host__

#endif

   void computeSpaceSteps();



   CoordinatesType dimensions;

src/mesh/tnlGrid2D_impl.h

+3 −0
Original line number Diff line number Diff line
@@ -76,6 +76,9 @@ tnlString tnlGrid< 2, Real, Device, Index > :: getSerializationTypeVirtual() con 
template< typename Real,

          typename Device,

          typename Index >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

void tnlGrid< 2, Real, Device, Index > :: computeSpaceSteps()

{

   if( this->getDimensions().x() > 0 && this->getDimensions().y() > 0 )

src/solvers/linear/krylov/tnlGMRESSolver_impl.h

+22 −2
Original line number Diff line number Diff line
@@ -105,7 +105,11 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
           << ". Please set some positive value using the SetRestarting method." << endl;

      return false;

   }

   if( ! setSize( matrix -> getRows(), restarting ) ) return false;

   if( ! setSize( matrix -> getRows(), restarting ) )

   {

       cerr << "I am not able to allocate enough memory for the GMRES solver. You may try to decrease the restarting parameter." << endl;

       return false;

   }





   IndexType i, j = 1, k, l;
@@ -145,6 +149,8 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
      normb = b. lpNorm( ( RealType ) 2.0 );

      _r. alphaXPlusBetaY( ( RealType ) 1.0, b, -1.0 );

      beta = _r. lpNorm( ( RealType ) 2.0 );

      cout << "x = " << x << endl;

      cout << " beta = " << beta << endl;

   }



   if( normb == 0.0 ) normb = 1.0;
@@ -195,6 +201,8 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
         else

             matrix -> vectorProduct( vi, _w );

         

         cout << " i = " << i << " vi = " << vi << endl;



         for( k = 0; k <= i; k++ )

         {

            vk. bind( &( _v. getData()[ k * _size ] ), _size );
@@ -208,6 +216,9 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
             * w = w - H_{k,i} v_k

             */

            _w. addVector( vk, -H_k_i );

            

            cout << "H_ki = " << H_k_i << endl;

            cout << "w = " << _w << endl;

         }

         /***

          * H_{i+1,i} = |w|
@@ -215,12 +226,16 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
         RealType normw = _w. lpNorm( ( RealType ) 2.0 );

         H[ i + 1 + i * ( m + 1 ) ] = normw;



         cout << "normw = " << normw << endl;

         

         /***

          * v_{i+1} = w / |w|

          */

         vi. bind( &( _v. getData()[ ( i + 1 ) * size ] ), size );

         vi. addVector( _w, ( RealType ) 1.0 / normw );

         

         cout << "vi = " << vi << endl;



         /****

          * Applying the Givens rotations

          */
@@ -276,6 +291,11 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
         beta = _r. lpNorm( ( RealType ) 2.0 );

      }

      this->setResidue( beta / normb );



      cout << " x = " << x << endl;

      cout << " beta = " << beta << endl;

      cout << "residue = " << beta / normb << endl;



   }

   this->refreshSolverMonitor();

   return this->checkConvergence();