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

Adding parameter tau to GetExplicitRHS in ODE solvers.

parent e517d66f

src/solvers/tnlEulerSolver.h

+1 −1
Original line number Diff line number Diff line
@@ -108,7 +108,7 @@ bool tnlEulerSolver< Problem, Mesh, Real, Device, Index > :: solve( Problem& sch 
      /****

       * Compute the RHS

       */

      scheme. GetExplicitRHS( time, u, k1 );

      scheme. GetExplicitRHS( time, currentTau, u, k1 );



      Real lastResidue = residue;

      computeNewTimeLevel( u, currentTau, residue );

src/solvers/tnlMersonSolver.h

+10 −10
Original line number Diff line number Diff line
@@ -300,35 +300,35 @@ void tnlMersonSolver< Problem, Mesh, Real, Device, Index > :: computeKFunctions( 


   if( Device :: getDevice() == tnlHostDevice )

   {

      problem. GetExplicitRHS( time, u, k1 );

      problem. GetExplicitRHS( time, tau, u, k1 );



   #ifdef HAVE_OPENMP

   #pragma omp parallel for firstprivate( size, _kAux, _u, _k1, tau, tau_3 )

   #endif

      for( Index i = 0; i < size; i ++ )

         _kAux[ i ] = _u[ i ] + tau * ( 1.0 / 3.0 * _k1[ i ] );

      problem. GetExplicitRHS( time + tau_3, kAux, k2 );

      problem. GetExplicitRHS( time + tau_3, tau, kAux, k2 );



   #ifdef HAVE_OPENMP

   #pragma omp parallel for firstprivate( size, _kAux, _u, _k1, _k2, tau, tau_3 )

   #endif

      for( Index i = 0; i < size; i ++ )

         _kAux[ i ] = _u[ i ] + tau * 1.0 / 6.0 * ( _k1[ i ] + _k2[ i ] );

      problem. GetExplicitRHS( time + tau_3, kAux, k3 );

      problem. GetExplicitRHS( time + tau_3, tau, kAux, k3 );



   #ifdef HAVE_OPENMP

   #pragma omp parallel for firstprivate( size, _kAux, _u, _k1, _k3, tau, tau_3 )

   #endif

      for( Index i = 0; i < size; i ++ )

         _kAux[ i ] = _u[ i ] + tau * ( 0.125 * _k1[ i ] + 0.375 * _k3[ i ] );

      problem. GetExplicitRHS( time + 0.5 * tau, kAux, k4 );

      problem. GetExplicitRHS( time + 0.5 * tau, tau, kAux, k4 );



   #ifdef HAVE_OPENMP

   #pragma omp parallel for firstprivate( size, _kAux, _u, _k1, _k3, _k4, tau, tau_3 )

   #endif

      for( Index i = 0; i < size; i ++ )

         _kAux[ i ] = _u[ i ] + tau * ( 0.5 * _k1[ i ] - 1.5 * _k3[ i ] + 2.0 * _k4[ i ] );

      problem. GetExplicitRHS( time + tau, kAux, k5 );

      problem. GetExplicitRHS( time + tau, tau, kAux, k5 );

   }

   if( Device :: getDevice() == tnlCudaDevice )

   {
@@ -336,27 +336,27 @@ void tnlMersonSolver< Problem, Mesh, Real, Device, Index > :: computeKFunctions( 
      const int block_size = 512;

      const int grid_size = ( size - 1 ) / block_size + 1;



      problem. GetExplicitRHS( time, u, k1 );

      problem. GetExplicitRHS( time, tau, u, k1 );

      cudaThreadSynchronize();



      computeK2Arg<<< grid_size, block_size >>>( size, tau, _u, _k1, _kAux );

      cudaThreadSynchronize();

      problem. GetExplicitRHS( time + tau_3, kAux, k2 );

      problem. GetExplicitRHS( time + tau_3, tau, kAux, k2 );

      cudaThreadSynchronize();



      computeK3Arg<<< grid_size, block_size >>>( size, tau, _u, _k1, _k2, _kAux );

      cudaThreadSynchronize();

      problem. GetExplicitRHS( time + tau_3, kAux, k3 );

      problem. GetExplicitRHS( time + tau_3, tau, kAux, k3 );

      cudaThreadSynchronize();



      computeK4Arg<<< grid_size, block_size >>>( size, tau, _u, _k1, _k3, _kAux );

      cudaThreadSynchronize();

      problem. GetExplicitRHS( time + 0.5 * tau, kAux, k4 );

      problem. GetExplicitRHS( time + 0.5 * tau, tau, kAux, k4 );

      cudaThreadSynchronize();



      computeK5Arg<<< grid_size, block_size >>>( size, tau, _u, _k1, _k3, _k4, _kAux );

      cudaThreadSynchronize();

      problem. GetExplicitRHS( time + tau, kAux, k5 );

      problem. GetExplicitRHS( time + tau, tau, kAux, k5 );

      cudaThreadSynchronize();

#endif

   }