TNL::Solvers::ODE::ODESolver< Method, Vector, SolverMonitor, IsStatic > Struct Template Reference

Integrator or solver of systems of ordinary differential equations. More...

Detailed Description

template<typename Method, typename Vector, typename SolverMonitor = IterativeSolverMonitor< typename Vector::RealType >, bool IsStatic = IsStaticArrayType< Vector >()>
struct TNL::Solvers::ODE::ODESolver< Method, Vector, SolverMonitor, IsStatic >

Integrator or solver of systems of ordinary differential equations.

This solver can be used for the numerical solution of ordinary differential equations having the following form:

\( \frac{d \vec u}{dt} = \vec f( t, \vec u) \text{ on } (0,T), \)

and the initial condition

\( \vec u( 0 ) = \vec u_{ini} \).

The vector \( \vec u(t) \) can be represented using different types of containers, depending on the size and nature of the ODE system:

1. Static vectors (TNL::Containers::StaticVector): This is suitable for small systems of ODEs with a fixed number of unknowns. Utilizing StaticVector allows the ODE solver to be executed within GPU kernels. This capability is particularly useful for scenarios like running multiple sequential solvers in parallel, as in the case of TNL::Algorithms::parallelFor.
2. Dynamic vectors (TNL::Containers::Vector or TNL::Containers::VectorView): These are preferred when dealing with large systems of ODEs, such as those arising in the solution of parabolic partial differential equations using the method of lines. In these instances, the solver typically handles a single, large-scale problem that can be executed in parallel internally.

The method, which is supposed to be used by the solver, is represented by the template parameter Method.

The following examples demonstrates the use the solver with the static vector

#include <iostream>
#include <fstream>
#include <TNL/Containers/Vector.h>
#include <TNL/Algorithms/parallelFor.h>
#include <TNL/Containers/StaticVector.h>
#include <TNL/Solvers/ODE/ODESolver.h>
#include <TNL/Solvers/ODE/Methods/Euler.h>

using MultiIndex = TNL::Containers::StaticArray< 3, int >;
using Real = double;
 
template< typename Device >
void
solveParallelODEs( const char* file_name )
{
   using StaticVector = TNL::Containers::StaticVector< 3, Real >;
   using Method = TNL::Solvers::ODE::Methods::Euler< Real >;
   using ODESolver = TNL::Solvers::ODE::ODESolver< Method, StaticVector >;

   const Real final_t = 50.0;
   const Real tau = 0.001;
   const Real output_time_step = 0.005;

   const Real sigma_min( 10.0 ), rho_min( 15.0 ), beta_min( 1.0 );
   const int parametric_steps = 4;
   const Real sigma_step = 30.0 / ( parametric_steps - 1 );
   const Real rho_step = 21.0 / ( parametric_steps - 1 );
   const Real beta_step = 15.0 / ( parametric_steps - 1 );
   const int output_time_steps = ceil( final_t / output_time_step ) + 1;

   const int results_size( output_time_steps * parametric_steps * parametric_steps * parametric_steps );
   TNL::Containers::Vector< StaticVector, Device > results( results_size, 0.0 );
   auto results_view = results.getView();

   auto f = [ = ] __cuda_callable__( const Real& t,
                                     const Real& tau,
                                     const StaticVector& u,
                                     StaticVector& fu,
                                     const Real& sigma_i,
                                     const Real& rho_j,
                                     const Real& beta_k )
   {
      const Real& x = u[ 0 ];
      const Real& y = u[ 1 ];
      const Real& z = u[ 2 ];
      fu[ 0 ] = sigma_i * ( y - x );
      fu[ 1 ] = rho_j * x - y - x * z;
      fu[ 2 ] = -beta_k * z + x * y;
   };

   auto solve = [ = ] __cuda_callable__( const MultiIndex& idx ) mutable
   {
      const Real sigma_i = sigma_min + idx[ 0 ] * sigma_step;
      const Real rho_j = rho_min + idx[ 1 ] * rho_step;
      const Real beta_k = beta_min + idx[ 2 ] * beta_step;

      ODESolver solver;
      solver.setTau( tau );
      solver.setTime( 0.0 );
      StaticVector u( 1.0, 1.0, 1.0 );
      int time_step( 1 );
      results_view[ ( idx[ 0 ] * parametric_steps + idx[ 1 ] ) * parametric_steps + idx[ 2 ] ] = u;

      while( time_step < output_time_steps ) {
         solver.setStopTime( TNL::min( solver.getTime() + output_time_step, final_t ) );
         solver.solve( u, f, sigma_i, rho_j, beta_k );
         const int k =
            ( ( time_step++ * parametric_steps + idx[ 0 ] ) * parametric_steps + idx[ 1 ] ) * parametric_steps + idx[ 2 ];
         results_view[ k ] = u;
      }
   };

   const MultiIndex begin = { 0, 0, 0 };
   const MultiIndex end = { parametric_steps, parametric_steps, parametric_steps };
   TNL::Algorithms::parallelFor< Device >( begin, end, solve );

   std::fstream file;
   file.open( file_name, std::ios::out );
   for( int sigma_idx = 0; sigma_idx < parametric_steps; sigma_idx++ )
      for( int rho_idx = 0; rho_idx < parametric_steps; rho_idx++ )
         for( int beta_idx = 0; beta_idx < parametric_steps; beta_idx++ ) {
            Real sigma = sigma_min + sigma_idx * sigma_step;
            Real rho = rho_min + rho_idx * rho_step;
            Real beta = beta_min + beta_idx * beta_step;
            file << "# sigma " << sigma << " rho " << rho << " beta " << beta << std::endl;
            for( int i = 0; i < output_time_steps - 1; i++ ) {
               int offset = ( ( i * parametric_steps + sigma_idx ) * parametric_steps + rho_idx ) * parametric_steps + beta_idx;
               auto u = results.getElement( offset );
               file << u[ 0 ] << " " << u[ 1 ] << " " << u[ 2 ] << std::endl;
            }
            file << std::endl;
         }
}
 
int
main( int argc, char* argv[] )
{
   if( argc != 2 ) {
      std::cout << "Usage: " << argv[ 0 ] << " <path to output directory>" << std::endl;
      return EXIT_FAILURE;
   }
   TNL::String file_name( argv[ 1 ] );
   file_name += "/StaticODESolver-LorenzParallelExample-result.out";
 
   solveParallelODEs< TNL::Devices::Host >( file_name.getString() );
#ifdef __CUDACC__
   solveParallelODEs< TNL::Devices::Cuda >( file_name.getString() );
#endif
}

and with the dynamic vector

#include <iostream>
#include <TNL/Containers/Vector.h>
#include <TNL/Solvers/ODE/ODESolver.h>
#include <TNL/Solvers/ODE/Methods/Euler.h>
#include "write.h"
 
using Real = double;
using Index = int;
 
template< typename Device >
void
solveHeatEquation( const char* file_name )
{
   using Vector = TNL::Containers::Vector< Real, Device, Index >;
   using VectorView = typename Vector::ViewType;
   using Method = TNL::Solvers::ODE::Methods::Euler< Real >;
   using ODESolver = TNL::Solvers::ODE::ODESolver< Method, Vector >;
 
   /***
    * Parameters of the discretisation
    */
   const Real final_t = 0.05;
   const Real output_time_step = 0.005;
   const Index n = 41;
   const Real h = 1.0 / ( n - 1 );
   const Real tau = 0.1 * h * h;
   const Real h_sqr_inv = 1.0 / ( h * h );
 
   /***
    * Initial condition
    */
   Vector u( n );
   u.forAllElements(
      [ = ] __cuda_callable__( Index i, Real & value )
      {
         const Real x = i * h;
         if( x >= 0.4 && x <= 0.6 )
            value = 1.0;
         else
            value = 0.0;
      } );
   std::fstream file;
   file.open( file_name, std::ios::out );
   write( file, u, n, h, (Real) 0.0 );
 
   /***
    * Setup of the solver
    */
   ODESolver solver;
   solver.setTau( tau );
   solver.setTime( 0.0 );
 
   /***
    * Time loop
    */
   while( solver.getTime() < final_t ) {
      solver.setStopTime( TNL::min( solver.getTime() + output_time_step, final_t ) );
      auto f = [ = ] __cuda_callable__( Index i, const VectorView& u, VectorView& fu ) mutable
      {
         if( i == 0 || i == n - 1 )  // boundary nodes -> boundary conditions
            fu[ i ] = 0.0;
         else  // interior nodes -> approximation of the second derivative
            fu[ i ] = h_sqr_inv * ( u[ i - 1 ] - 2.0 * u[ i ] + u[ i + 1 ] );
      };
      auto time_stepping = [ = ]( const Real& t, const Real& tau, const VectorView& u, VectorView& fu )
      {
         TNL::Algorithms::parallelFor< Device >( 0, n, f, u, fu );
      };
      solver.solve( u, time_stepping );
      write( file, u, n, h, solver.getTime() );  // write the current state to a file
   }
}
 
int
main( int argc, char* argv[] )
{
   if( argc != 2 ) {
      std::cout << "Usage: " << argv[ 0 ] << " <path to output directory>" << std::endl;
      return EXIT_FAILURE;
   }
   TNL::String file_name( argv[ 1 ] );
   file_name += "/ODESolver-HeatEquationExample-result.out";
 
   solveHeatEquation< TNL::Devices::Host >( file_name.getString() );
#ifdef __CUDACC__
   solveHeatEquation< TNL::Devices::Cuda >( file_name.getString() );
#endif
}

Template Parameters

Method	is a method (one from TNL::Solvers::ODE namespace) which is supposed to be used for the numerical integration.
Vector	is a vector (TNL::Containers::Vector, TNL::Containers::VectorView, or TNL::Containers::StaticVector) representing \( \vec x \in R^n \).

---

The documentation for this struct was generated from the following file:

• src/TNL/Solvers/ODE/ODESolver.h