Table of Contents
Introduction
In this part, we describe solvers of ordinary differential equations having the following form:
\[ \frac{d \vec u}{dt} = \vec f( t, \vec u) \text{ on } (0,T), \]
\[ \vec u( 0 ) = \vec u_{ini}. \]
TNL offers the following ODE solvers:
- TNL::Solvers::ODE::Euler - the Euler method with the 1-st order of accuracy.
- TNL::Solvers::ODE::Merson - the Runge-Kutta-Merson solver with the 4-th order of accuracy and adaptive choice of the time step.
Each solver has its static counterpart which can be run even in the GPU kernels which means that it can be combined with TNL::Algorithms::parallelFor for example. The static ODE solvers are the following:
- TNL::Solvers::ODE::StaticEuler - the Euler method with the 1-st order of accuracy.
- TNL::Solvers::ODE::StaticMerson - the Runge-Kutta-Merson solver with the 4-th order of accuracy and adaptive choice of the time step.
Static ODE solvers
Static solvers are supposed to be used mainly when \( x \in R \) is scalar or \( \vec x \in R^n \) is vector where \( n \) is small. Firstly, we show example of scalar problem of the following form:
\[ \frac{d u}{dt} = t \sin ( t )\ \rm{ on }\ (0,T), \]
\[ u( 0 ) = 0. \]
This problem can be solved as follows:
We first define the type Real representing the floating-point arithmetics, which is double in this case (line 4). In the main function, we define the type of the ODE solver (ODESolver, line 8). We choose TNL::Solvers::ODE::StaticEuler. We define the variable final_t (line 9) representing the size of the time interval \( (0,T)\), next we define the integration time step tau (line 10) and the variable output_time_step (line 11) representing checkpoints in which we will print value of the solution \( x(t)\). On the line 13, we create an instance of the ODESolver and set the integration time step (line 14) and the initial time of the solver (line 15). Next we create variable u representing the solution of the given ODE and we initiate it with the initial condition \( u(0) = 0\) (line 16). On the lines 17-25, we iterate over the interval \( (0,T) \) with step given by the frequency of the checkpoints given by the variable output_time_steps. On the line 19, we let the solver to iterate until the next checkpoint or the end of the interval \((0,T) \) depending on what occurs first (it is expressed by TNL::min( solver.getTime() + output_time_step, final_t )). On the lines 20-22, we define the lambda function f representing the right-hand side of the ODE being solved. The lambda function receives the following arguments:
tis the current value of the time variable \( t \in (0,T)\),tauis the current integration time step,uis the current value of the solution \( u(t)\),fuis a reference on a variable into which we evaluate the right-hand side \( f(u,t) \) on the ODE.
The lambda function is supposed to compute just the value of fu. It is t * sin(t) on our case. Finally we call the ODE solver (line 23). As parameters, we pass the variable u representing the solution \( u(t)\) and a lambda function representing the right-hand side of the ODE. On the line 23, we print values of the solution at given checkpoints. The result looks as follows:
Such data can by visualized using Gnuplot as follows
or it can be processed by the following Python script which draws graph of the function \( u(t) \) using Matplotlib.
We first parse the input file (lines 13-20) and convert the data into NumPy arrays (lines 24-25). Finaly, these arrays are drawn into a graph (lines 29-36). The graph of the solution looks as follows:
In the next example, we demonstrate use of the static ODE solver to solve a system of ODEs, namely the Lorenz system which reads as:
\[ \frac{dx}{dt} = \sigma( x - y),\ \rm{ on }\ (0,T) \]
\[ \frac{dy}{dt} = x(\rho - z ) - y,\ \rm{ on }\ (0,T) \]
\[ \frac{dz}{dt} = xy - \beta z,\ \rm{ on }\ (0,T) \]
\[ \vec u(0) = (x(0),y(0),z(0)) = \vec u_{ini} \]
for given constants \( \sigma, \rho \) and \( \beta \). The solution \( \vec u(t) = (x(t), y(t), z(t)) \in R^3 \) is represented by three-dimensional static vector (TNL::Containers::StaticVector). The solver looks as follows:
The code is very similar to the previous example. There are the following differences:
- We define the type of the variable
urepresenting the solution \( \vec u(t) \) as TNL::Containers::StaticVector< 3, Real > (line 9) which is reflected even in the definition of the ODE solver (TNL::Solvers::ODE::StaticEuler< Vector > , line 10) and the variableu(line 21). - In addition to the parameters of the solver (
final_t,tauandoutput_time_step, lines 14-16) we define parameters of the Lorenz system (sigma,rhoandbeta, lines 14-16). - The initial condition \( \vec u(0) = (1,2,3) \) is set on the line 21.
- In the lambda function representing the right-hand side of the Lorenz system (lines 25-32), we first define auxiliary aliases
x,yandz(lines 26-28) to make the code easier to read. The main right-hand side of the Lorenz system is implemented on the lines (29-31). - In the line 3, we print all components of the vector
u.
The solver creates file with the solution \( (\sigma(i \tau), \rho( i \tau), \beta( i \tau )) \) for ' \( i = 0, 1, \ldots N \) on separate lines. It looks as follows:
Such file can visualized using Gnuplot interactively in 3D as follows
or it can be processed by the following Python script:
The script has very similar structure as in the previous example. The result looks as follows:
Combining static ODE solvers with parallel for
The static solvers can be used inside of lambda functions for TNL::Algorithms::parallelFor for example. This can be useful when we need to solve large number of independent ODE problems, for example for parametric analysis. We demonstrate it on the two examples we have described above.
Solving scalar problems in parallel
The first example solves ODE given by
\[ \frac{d u}{dt} = t \sin ( c t )\ \rm{ on }\ (0,T), \]
\[ u( 0 ) = 0, \]
where \( c \) is a constant. We will solve it in parallel ODEs with different values \( c \in \langle c_{min}, c_{max} \rangle \). The exact solution can be found here. The code reads as follows:
In this example we also show, how to run it on GPU. Therefore we moved the main solver to separate function solveParallelODEs which has one template parameter Device telling on what device it is supposed to run. The results of particular ODEs are stored in a memory and at the end they are copied to a file with given filename file_name. The variable \( u \) is scalar therefore we represent it by the type Real in the solver (line 12). Next we define parameters of the ODE solver (final_t, tau and output_time_step, lines 14-16), interval for the parameter of the ODE \( c \in \langle c_{min}, c_{max} \rangle \) ( c_min, c_max, lines 17-18) and number of values c_vals (line 19) distributed equidistantly in the interval with step c_step (line 20). We use the number of different values of the parameter c as a range for the parallel for on the line 43. This parallel for processes the lambda function solve which is defined on the lines 28-42. It receives a parameter idx which is index of the value of the parameter c. We compute its value on the line 29. Next we create the ODE solver (line 30) and setup its parameters (lines 31-32). We set initial condition of the given ODE and we define variable time_step which counts checkpoints which we store in the memory in vector results allocated in the line 23 and accessed in the lambda function via vector view results_view (defined on the line 24). We iterate over the interval \( (0, T) \) in the while loop starting on the line 36. We set the stop time of the ODE solver (line 38) and we run the solver (line 39). Finally we store the result at given checkpoint into vector view results_view. If the solver runs on the GPU, it cannot write the checkpoints into a file. This is done in postprocessing in the lines 45-53.
Note, how we pass the value of parameter c to the lambda function f. The method solve of the ODE solvers(TNL::Solvers::ODE::StaticEuler::solve, for example) accepts user defined parameters via variadic templates. It means that in addition to the variable u and the right-hand side f we can add any other parameters like c in this example (line 39). This parameter appears in the lambda function f (line 25). The reason for this is that the nvcc compiler (version 10.1) does not accept lambda function defined within another lambda function. If such a construction is accepted by a compiler, the lambda function f which can be defined within the lambda function solve and the variable c defined in the lambda function solve could be captured by f.
The solver generates file of the following format:
The file an visualized using Gnuplot as follows
or it can be processed by the following Python script:
The result of this example looks as follows:
Solving the Lorenz system in parallel
The second example is a parametric analysis of the Lorenz model
\[ \frac{dx}{dt} = \sigma( x - y),\ \rm{ on }\ (0,T) \]
\[ \frac{dy}{dt} = x(\rho - z ) - y,\ \rm{ on }\ (0,T) \]
\[ \frac{dz}{dt} = xy - \beta z,\ \rm{ on }\ (0,T) \]
\[ \vec u(0) = (x(0),y(0),z(0)) = \vec u_{ini} \]
which we solve for different values of the model parameters
\[ \sigma_i = \sigma_{min} + i \Delta \sigma, \]
\[ \rho_j = \rho_{min} + j \Delta \rho, \]
\[ \beta_k = \beta_{min} + k \Delta \beta, \]
where we set \( \Delta \sigma = \Delta \rho = \Delta \beta = l / (p-1) \) and where \( i,j,k = 0, 1, \ldots, p - 1 \). The code of the solver looks as follows:
It is very similar to the previous one. There are just the following changes:
- On the line 17, we define minimal values for the parameters \( \sigma, \beta \) and \( \rho \). On the line 18, we define how many different values we will consider for each parameter. The size of equidistant steps in the parameter variations is defined on the line 19. The interval for parameters variations is set to \( [0,30] \) (line 19 as well).
- On the line 23, we allocate vector
resultsinto which we will store the solution of the Lorenz problem for various parameters. - Next we define the lambda function
frepresenting the right-hand side of the Lorenz problem (lines 25-33) and the lambda functionsolverepresenting the ODE solver for the Lorenz problem with particular setup of the parameters (lines 34-52). This lambda function is processed by TNL::Algorithms::parallelFor called on the line 53. Therefore the lambda functionsolvereceives three indexesi,jandkwhich are used to compute particular values of the parameters \( \sigma_i, \rho_j, \beta_k \) which are represented by variablessigma_i,rho_jandbeta_k(lines 35-37). These parameters must be passed to the lambda functionfexplicitly (line 48). The reason is the same as in the previous example - nvcc (version 10.1) does not accept a lambda function defined within another lambda function. - The initial condition for the Lorenz problem is set to vector \( (1,1,1) \) (line 42). Finally, we start the time loop (lines 45-51) and we store the state of the solution into the vector
resultsusing related vector viewresults_viewin the time intervals given by the variableoutput_time_step. - When all ODEs are solved, we copy all the solutions from the vector
resultsinto an output file.
The files has the following format:
The file can be processed by the following Python script:
The result looks as follows:
Non-static ODE Solvers
In this section, we will show how to solve simple heat equation in 1D. The heat equation is a parabolic partial differential equation which reads as follows:
\[ \frac{\partial u(t,x)}{\partial t} - \frac{\partial^2 u(t,x)}{\partial^2 x} = 0\ \rm{on}\ (0,T) \times (0,1), \]
\[ u(t,0) = 0, \]
\[ u(t,0) = 1, \]
\[ u(0,x) = u_{ini}(x)\ \rm{on}\ [0,1]. \]
We discretize it by the finite difference method to get numerical approximation. We first define set of nodes \( x_i = ih \) for \(i=0,\ldots n-1 \) where \(h = 1 / (n-1) \) (we use C++ indexing for the consistency). Using the method of lines and approximating the second derivative by the central finite diference ( \( \frac{\partial^2 u(t,x)}{\partial^2 x} \approx \frac{u_{i-1} - 2 u_i + u_{i+1}}{h^2} \)), we obtain system of ODEs of the following form
\[ \frac{\rm{d}u_i(t)}{\rm{d} t} = \frac{u_{i-1} - 2 u_i + u_{i+1}}{h^2}\ \rm{for}\ i = 1, \ldots, n-2, \]
where \( u_i(t) = u(t,ih) \) and \( h \) space step between two adjacent nodes of the numerical mesh. We also set
\[ u_0(t) = u_{n-1}(t) = 0\ \rm{on}\ [0,T]. \]
What are the main differences compared to the Lorenz model?
- The size of the Lorenz model is fixed. It is equal to three since \( (\sigma, \rho, \beta) \in R^3 \) which is small vector of fixed size and it can be represented by the static vector TNL::Containers::StaticVector< 3, Real >. On the other hand, the size of the ODE system arising in the discretization by the method of lines depends not on the problem we solve but on the desired accuracy - the larger \( n \) the more accurate numerical approximation we get. The number of nodes \( n \) used for the space discretisation defines the number of parameters defining the mesh function. These parameters are also referred as degrees of freedom, DOFs. Therefore the size of the system can be large and it is better to employ dynamic vector TNL::Containers::Vector for the solution.
- The size of the Lorenz model is small and so the evaluation of its right-hand side can be done sequentially by one thread. The size of the ODE system can be very large and evaluating the right-hand side asks for being performed in parallel.
- The dynamic vector TNL::Containers::Vector allocates data dynamically and therefore it cannot be created within a GPU kernel which means the ODE solvers cannot be created in the GPU kernel either. For this reason, the lambda function
fevaluating the right-hand side of the ODE system is always executed on the host and it calls TNL::Algorithms::parallelFor to evaluate the right-hand side of the ODE system.
Basic setup
The implementation of the solver looks as follows:
The solver is implemented in separate function solveHeatEquation (line 11) having one template parameter Device (line 10) defining on what device the solver is supposed to be executed. We first define necessary types:
Vector(line 13) is alias for TNL::Containers::Vector. The number of DOFs can be large and therefore they are stored in the resizable dynamically allocated vector TNL::Containers::Vector rather then in static vector - TNL::Containers::StaticVector.VectorView(line 14) will be used for accessing DOFs within the lambda functions.ODESolver(line 15) is type of the ODE solver which we will use for the computation of the time evolution in the time dependent heat equation.
Next we define parameters of the discretization:
final_t(line 20) represents the size of the time interval \( [0,T] \).output_time_step(line 21) defines time intervals in which we will write the solution \( u \) into a file.n(line 22) is the number of DOFs, i.e. number of nodes used for the finite difference method.h(line 23) is the space step, i.e. distance between two consecutive nodes.tau(line 24) is the time step used for the integration by the ODE solver. Since we solve the second order parabolic problem, the time step should be proportional to \( h^2 \).h_sqr_inv(line 25) is auxiliary constant equal to \( 1/h^2 \) which we use later in finite difference for approximation of the second derivative.
Next we set the initial condition \( u_{ini} \) (lines ) which is given as:
\[ u_{ini}(x) = \left\{ \begin{array}{rl} 0 & \rm{for}\ x < 0.4, \\ 1 & \rm{for}\ 0.4 \leq x \leq 0.6, \\ 0 & \rm{for}\ x > 0. \\ \end{array} \right. \]
Next we write the initial condition to a file (lines 37-39) using the function write which we will describe later. On the lines (44-46) we create an instance of the ODE solver solver, we set the integration time step tau of the solver (TNL::Solvers::ODE::ExplicitSolver::setTau ) and we set the initial time to zero (TNL::Solvers::ODE::ExplicitSolver::setTime).
Finally, we proceed to the time loop (lines 51-64) but before we prepare counter of the states to be written into files (output_idx). The time loop uses the time variable within the ODE solver (TNL::Solvers::ODE::ExplicitSolver::getTime ) and it iterates until we reach the end of the time interval \( [0, T] \) given by the variable final_t. On the line 53, we set the stop time of the ODE solver (TNL::Solvers::ODE::ExplicitSolver::setStopTime ) to the next checkpoint for storing the state of the heat equation or the end of the time interval depending on what comes first. Next we define the lambda function f expressing the discretization of the second derivative of \( u \) by the central finite difference and the boundary conditions. The function receives the following parameters:
iis the index of the node and the related ODE arising from the method of lines. In fact, we have to evaluate the update of \( u_i^k \) to get to the next time level \( u_i^{k+1} \).uis vector view representing the state \( u_i^k \) of the heat equation on the \( k- \) time level.fuis vector of updates or time derivatives in the method of lines which will bring \( u \) to the next time level.
As we mentioned above, since nvcc does not accept lambda functions defined within another lambda function, we have to define f separately and pass the parameters u and fu explicitly (see the line 61).
Now look at the code of the lambda function f. Since the solution \( u \) does not change on the boundaries, we return zero on the boundary nodes (lines 55-56) and we evaluate the central difference for approximation of the second derivative on the interior nodes (line 58).
Next we define the lambda function time_stepping (lines ) which is responsible for computing of the updates for all nodes \( i = 0, \ldots n-1 \). It is done by means of TNL::Algorithms::parallelFor which iterates over all the nodes and calling the function f on each of them. It passes the vector views u and fu explicitly to f for the reasons we have mentioned above.
Finally, we run the ODE solver (TNL::Solvers::ODE::Euler::solve ) (line 62) and we pass u as the current state of the heat equation and f the lambda function controlling the time evolution to the method solve. On the line 63, we store the current state to a file.
The function write which we use for writing the solution of the heat equation reads as follows:
The function accepts the following parameters:
fileis the file into which we want to store the solution.uis a vector or vector view representing solution of the heat equation at given time.nis number of nodes the we use for the approximation of the solution.his space step, i.e. distance between two consecutive nodes.timeis the current time of the evolution being computed.
The solver writes the results in the following format:
The solution can be visualised with Gnuplot as follows:
or it can be easily parsed in Python and processes by Matplotlib using the following script:
The result looks as follows:
Setup with a solver monitor
In this section we will show how to connect ODE solver with the solver monitor. Look at the following example:
There are the following differences compared to the previous example:
- We have to include a header file with the iterative solver monitor (line 5).
- We have to setup the solver monitor (lines 45-50). First, we define the monitor type (line 45) and we create an instance of the monitor (line 46). Next we create a separate thread for the monitor (line 47), set the refresh rate to 10 milliseconds (line 48), turn on the verbose mode (line 49) and set the solver stage name (line 50). On the line 58, we connect the monitor with the solver using method TNL::Solvers::IterativeSolver::setSolverMonitor. We stop the monitor after the ODE solver finishes by calling TNL::Solvers::IterativeSolverMonitor::stopMainLoop in the line 77.
The result looks as follows:
