Template Numerical Library version\ main:978fbd2
Loading...
Searching...
No Matches
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, typename Vector::IndexType >, 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

1#include <iostream>
2#include <fstream>
3#include <TNL/Containers/Vector.h>
4#include <TNL/Algorithms/parallelFor.h>
5#include <TNL/Containers/StaticVector.h>
6#include <TNL/Solvers/ODE/ODESolver.h>
7#include <TNL/Solvers/ODE/Methods/Euler.h>
8
10using MultiIndex = TNL::Containers::StaticArray< 3, int >;
11using Real = double;
13
14template< typename Device >
15void
16solveParallelODEs( const char* file_name )
17{
19 using StaticVector = TNL::Containers::StaticVector< 3, Real >;
20 using Method = TNL::Solvers::ODE::Methods::Euler< Real >;
21 using ODESolver = TNL::Solvers::ODE::ODESolver< Method, StaticVector >;
23
25 const Real final_t = 50.0;
26 const Real tau = 0.001;
27 const Real output_time_step = 0.005;
29
31 const Real sigma_min( 10.0 ), rho_min( 15.0 ), beta_min( 1.0 );
32 const int parametric_steps = 4;
33 const Real sigma_step = 30.0 / ( parametric_steps - 1 );
34 const Real rho_step = 21.0 / ( parametric_steps - 1 );
35 const Real beta_step = 15.0 / ( parametric_steps - 1 );
36 const int output_time_steps = ceil( final_t / output_time_step ) + 1;
38
40 const int results_size( output_time_steps * parametric_steps * parametric_steps * parametric_steps );
41 TNL::Containers::Vector< StaticVector, Device > results( results_size, 0.0 );
42 auto results_view = results.getView();
44
46 auto f = [ = ] __cuda_callable__( const Real& t,
47 const Real& tau,
48 const StaticVector& u,
49 StaticVector& fu,
50 const Real& sigma_i,
51 const Real& rho_j,
52 const Real& beta_k )
53 {
54 const Real& x = u[ 0 ];
55 const Real& y = u[ 1 ];
56 const Real& z = u[ 2 ];
57 fu[ 0 ] = sigma_i * ( y - x );
58 fu[ 1 ] = rho_j * x - y - x * z;
59 fu[ 2 ] = -beta_k * z + x * y;
60 };
62
64 auto solve = [ = ] __cuda_callable__( const MultiIndex& idx ) mutable
65 {
67 const Real sigma_i = sigma_min + idx[ 0 ] * sigma_step;
68 const Real rho_j = rho_min + idx[ 1 ] * rho_step;
69 const Real beta_k = beta_min + idx[ 2 ] * beta_step;
71
73 ODESolver solver;
74 solver.setTau( tau );
75 solver.setTime( 0.0 );
76 StaticVector u( 1.0, 1.0, 1.0 );
77 int time_step( 1 );
78 results_view[ ( idx[ 0 ] * parametric_steps + idx[ 1 ] ) * parametric_steps + idx[ 2 ] ] = u;
80
82 while( time_step < output_time_steps ) {
83 solver.setStopTime( TNL::min( solver.getTime() + output_time_step, final_t ) );
85 solver.solve( u, f, sigma_i, rho_j, beta_k );
87 const int k =
88 ( ( time_step++ * parametric_steps + idx[ 0 ] ) * parametric_steps + idx[ 1 ] ) * parametric_steps + idx[ 2 ];
89 results_view[ k ] = u;
90 }
92 };
94
96 const MultiIndex begin = { 0, 0, 0 };
97 const MultiIndex end = { parametric_steps, parametric_steps, parametric_steps };
98 TNL::Algorithms::parallelFor< Device >( begin, end, solve );
100
102 std::fstream file;
103 file.open( file_name, std::ios::out );
104 for( int sigma_idx = 0; sigma_idx < parametric_steps; sigma_idx++ )
105 for( int rho_idx = 0; rho_idx < parametric_steps; rho_idx++ )
106 for( int beta_idx = 0; beta_idx < parametric_steps; beta_idx++ ) {
107 Real sigma = sigma_min + sigma_idx * sigma_step;
108 Real rho = rho_min + rho_idx * rho_step;
109 Real beta = beta_min + beta_idx * beta_step;
110 file << "# sigma " << sigma << " rho " << rho << " beta " << beta << std::endl;
111 for( int i = 0; i < output_time_steps - 1; i++ ) {
112 int offset = ( ( i * parametric_steps + sigma_idx ) * parametric_steps + rho_idx ) * parametric_steps + beta_idx;
113 auto u = results.getElement( offset );
114 file << u[ 0 ] << " " << u[ 1 ] << " " << u[ 2 ] << std::endl;
115 }
116 file << std::endl;
117 }
119}
120
121int
122main( int argc, char* argv[] )
123{
124 if( argc != 2 ) {
125 std::cout << "Usage: " << argv[ 0 ] << " <path to output directory>" << std::endl;
126 return EXIT_FAILURE;
127 }
128 TNL::String file_name( argv[ 1 ] );
129 file_name += "/StaticODESolver-LorenzParallelExample-result.out";
130
131 solveParallelODEs< TNL::Devices::Host >( file_name.getString() );
132#ifdef __CUDACC__
133 solveParallelODEs< TNL::Devices::Cuda >( file_name.getString() );
134#endif
135}
__cuda_callable__
#define __cuda_callable__
Definition Macros.h:49
TNL::Containers::StaticArray
Array with constant size.
Definition StaticArray.h:20
TNL::Containers::StaticVector
Vector with constant size.
Definition StaticVector.h:19
TNL::Containers::Vector
Vector extends Array with algebraic operations.
Definition Vector.h:36
TNL::Real
Definition Real.h:14
TNL::String
Class for managing strings.
Definition String.h:30
std::endl
T endl(T... args)
TNL::Algorithms::parallelFor
std::enable_if_t< std::is_integral_v< Begin > &&std::is_integral_v< End > > parallelFor(const Begin &begin, const End &end, typename Device::LaunchConfiguration launch_config, Function f, FunctionArgs... args)
Parallel for-loop function for 1D range specified with integral values.
Definition parallelFor.h:41
std::fstream::open
T open(T... args)
TNL::Solvers::ODE::Methods::Euler
First order Euler method.
Definition Euler.h:15
TNL::Solvers::ODE::ODESolver
Integrator or solver of systems of ordinary differential equations.
Definition ODESolver.h:61

and with the dynamic vector

1#include <iostream>
2#include <TNL/FileName.h>
3#include <TNL/Containers/Vector.h>
4#include <TNL/Solvers/ODE/ODESolver.h>
5#include <TNL/Solvers/ODE/Methods/Euler.h>
6#include "write.h"
7
8using Real = double;
9using Index = int;
10
11template< typename Device >
12void
13solveHeatEquation( const char* file_name )
14{
16 using Vector = TNL::Containers::Vector< Real, Device, Index >;
17 using VectorView = typename Vector::ViewType;
18 using Method = TNL::Solvers::ODE::Methods::Euler< Real >;
19 using ODESolver = TNL::Solvers::ODE::ODESolver< Method, Vector >;
21
22 /***
23 * Parameters of the discretisation
24 */
26 const Real final_t = 0.05;
27 const Real output_time_step = 0.005;
28 const Index n = 41;
29 const Real h = 1.0 / ( n - 1 );
30 const Real tau = 0.1 * h * h;
31 const Real h_sqr_inv = 1.0 / ( h * h );
33
34 /***
35 * Initial condition
36 */
38 Vector u( n );
39 u.forAllElements(
40 [ = ] __cuda_callable__( Index i, Real & value )
41 {
42 const Real x = i * h;
43 if( x >= 0.4 && x <= 0.6 )
44 value = 1.0;
45 else
46 value = 0.0;
47 } );
48 std::fstream file;
49 file.open( file_name, std::ios::out );
50 write( file, u, n, h, (Real) 0.0 );
52
53 /***
54 * Setup of the solver
55 */
57 ODESolver solver;
58 solver.setTau( tau );
59 solver.setTime( 0.0 );
61
62 /***
63 * Time loop
64 */
66 while( solver.getTime() < final_t ) {
67 solver.setStopTime( TNL::min( solver.getTime() + output_time_step, final_t ) );
69 auto f = [ = ] __cuda_callable__( Index i, const VectorView& u, VectorView& fu ) mutable
70 {
71 if( i == 0 || i == n - 1 ) // boundary nodes -> boundary conditions
72 fu[ i ] = 0.0;
73 else // interior nodes -> approximation of the second derivative
74 fu[ i ] = h_sqr_inv * ( u[ i - 1 ] - 2.0 * u[ i ] + u[ i + 1 ] );
75 };
78 auto time_stepping = [ = ]( const Real& t, const Real& tau, const VectorView& u, VectorView& fu )
79 {
81 TNL::Algorithms::parallelFor< Device >( 0, n, f, u, fu );
83 };
85 solver.solve( u, time_stepping );
86 write( file, u, n, h, solver.getTime() ); // write the current state to a file
87 }
89}
90
91int
92main( int argc, char* argv[] )
93{
94 if( argc != 2 ) {
95 std::cout << "Usage: " << argv[ 0 ] << " <path to output directory>" << std::endl;
96 return EXIT_FAILURE;
97 }
98 TNL::String file_name( argv[ 1 ] );
99 file_name += "/ODESolver-HeatEquationExample-result.out";
100
101 solveHeatEquation< TNL::Devices::Host >( file_name.getString() );
102#ifdef __CUDACC__
103 solveHeatEquation< TNL::Devices::Cuda >( file_name.getString() );
104#endif
105}
Template Parameters
Methodis a method (one from TNL::Solvers::ODE namespace) which is supposed to be used for the numerical integration.
Vectoris 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:
Generated on Tue Aug 27 2024 13:11:49 for Template Numerical Library by doxygen 1.12.0