Skip to content
Snippets Groups Projects
Commit c411380f authored by Tomáš Oberhuber's avatar Tomáš Oberhuber Committed by Jakub Klinkovský
Browse files

Added Python scripts for static ODE solver examples.

parent a19f4181
No related branches found
No related tags found
1 merge request!125ODE solvers
......@@ -11,22 +11,30 @@ set( COMMON_EXAMPLES
)
foreach( target IN ITEMS ${HOST_EXAMPLES})
add_executable( ${target} ${target}.cpp )
add_custom_command( COMMAND ${target} > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out OUTPUT ${target}.out )
set( HOST_OUTPUTS ${HOST_OUTPUTS} ${target}.out )
add_executable( ${target} ${target}.cpp )
add_custom_command( COMMAND ${target} > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out
OUTPUT ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}-result.out )
set( HOST_OUTPUTS ${HOST_OUTPUTS} ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out )
endforeach()
if( BUILD_CUDA )
foreach( target IN ITEMS ${COMMON_EXAMPLES} )
cuda_add_executable( ${target}-cuda ${target}.cu OPTIONS )
add_custom_command( COMMAND ${target}-cuda > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out OUTPUT ${target}.out )
set( CUDA_OUTPUTS ${CUDA_OUTPUTS} ${target}.out )
add_custom_command( COMMAND ${target}-cuda ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}
> ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out
OUTPUT ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}-result.out)
set( CUDA_OUTPUTS ${CUDA_OUTPUTS} ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out )
endforeach()
else()
foreach( target IN ITEMS ${COMMON_EXAMPLES})
foreach( target IN ITEMS ${COMMON_EXAMPLES} )
add_executable( ${target} ${target}.cpp )
add_custom_command( COMMAND ${target} > ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out OUTPUT ${target}.out )
set( HOST_OUTPUTS ${HOST_OUTPUTS} ${target}.out )
add_custom_command( COMMAND ${target} ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}
> ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out
OUTPUT ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}-result.out )
set( HOST_OUTPUTS ${HOST_OUTPUTS} ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/${target}.out )
endforeach()
endif()
......@@ -44,6 +52,14 @@ ADD_CUSTOM_COMMAND(
DEPENDS ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-SineExample.out
VERBATIM )
ADD_CUSTOM_COMMAND(
COMMAND python3 ${CMAKE_CURRENT_SOURCE_DIR}/StaticODESolver-SineParallelExample.py
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-SineParallelExample-result.out
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-SineParallelExample.png
OUTPUT StaticODESolver-SineParallelExample.png
DEPENDS ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-SineParallelExample.out
VERBATIM )
ADD_CUSTOM_COMMAND(
COMMAND python3 ${CMAKE_CURRENT_SOURCE_DIR}/StaticODESolver-LorenzExample.py
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-LorenzExample.out
......@@ -52,6 +68,17 @@ ADD_CUSTOM_COMMAND(
DEPENDS ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-LorenzExample.out
VERBATIM )
ADD_CUSTOM_TARGET( ProcessODESolversExamplesResults ALL DEPENDS
StaticODESolver-SineExample.png
StaticODESolver-LorenzExample.png )
\ No newline at end of file
ADD_CUSTOM_COMMAND(
COMMAND python3 ${CMAKE_CURRENT_SOURCE_DIR}/StaticODESolver-LorenzParallelExample.py
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-LorenzParallelExample-result.out
${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-LorenzParallelExample
OUTPUT StaticODESolver-LorenzParallelExample.png
DEPENDS ${TNL_DOCUMENTATION_OUTPUT_SNIPPETS_PATH}/StaticODESolver-LorenzParallelExample.out
VERBATIM )
ADD_CUSTOM_TARGET( ProcessODESolversExamplesResults ALL DEPENDS
StaticODESolver-SineExample.png
StaticODESolver-SineParallelExample.png
StaticODESolver-LorenzExample.png
StaticODESolver-LorenzParallelExample.png )
\ No newline at end of file
......@@ -2,6 +2,7 @@
import sys
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D
import numpy as np
plt.rcParams['text.usetex'] = True
......@@ -22,7 +23,8 @@ y = np.array(y_lst)
z = np.array(z_lst)
fig = plt.figure()
ax = fig.gca(projection='3d')
#ax = fig.gca(projection='3d')
ax = Axes3D(fig)
theta = np.linspace(-4 * np.pi, 4 * np.pi, 100)
ax.plot(x, y, z, label='Lorenz attractor')
ax.legend()
......
......@@ -11,12 +11,14 @@ void solveParallelODEs( const char* file_name )
{
using Vector = TNL::Containers::StaticVector< 3, Real >;
using ODESolver = TNL::Solvers::ODE::StaticEuler< Vector >;
const Real final_t = 10.0;
const Real final_t = 50.0;
const Real tau = 0.001;
const Real output_time_step = 0.01;
const Real sigma_min( 0.0 ), rho_min( 0.0 ), beta_min( 0.0 );
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 parametric_step = 30.0 / ( parametric_steps - 1 );
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 = final_t / output_time_step + 1;
const int results_size( output_time_steps * parametric_steps * parametric_steps * parametric_steps );
......@@ -32,16 +34,16 @@ void solveParallelODEs( const char* file_name )
fu[ 2 ] = -beta_k * z + x * y;
};
auto solve = [=] ( int i, int j, int k ) mutable {
const Real sigma_i = sigma_min + i * parametric_step;
const Real rho_j = rho_min + j * parametric_step;
const Real beta_k = beta_min + k * parametric_step;
const Real sigma_i = sigma_min + i * sigma_step;
const Real rho_j = rho_min + j * rho_step;
const Real beta_k = beta_min + k * beta_step;
ODESolver solver;
solver.setTau( tau );
solver.setTime( 0.0 );
Vector u( 1.0, 1.0, 1.0 );
int time_step( 0 );
results_view[ ( i * parametric_step + j ) * parametric_step + k ] = u;
results_view[ ( i * parametric_steps + j ) * parametric_steps + k ] = u;
while( solver.getTime() < final_t )
{
solver.setStopTime( TNL::min( solver.getTime() + output_time_step, final_t ) );
......@@ -58,11 +60,11 @@ void solveParallelODEs( const char* file_name )
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 * parametric_step;
Real rho = rho_min + rho_idx * parametric_step;
Real beta = beta_min + beta_idx * parametric_step;
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; i++ )
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 );
......@@ -74,8 +76,11 @@ void solveParallelODEs( const char* file_name )
int main( int argc, char* argv[] )
{
solveParallelODEs< TNL::Devices::Host >( "StaticODESolver-LorenzParallelExample-Host.out" );
TNL::String file_name( argv[ 1 ] );
file_name += "/StaticODESolver-LorenzParallelExample-result.out";
solveParallelODEs< TNL::Devices::Host >( file_name.getString() );
#ifdef HAVE_CUDA
solveParallelODEs< TNL::Devices::Cuda >( "StaticODESolver-LorenzParallelExample-Cuda.out" );
solveParallelODEs< TNL::Devices::Cuda >( file_name.getString() );
#endif
}
#!/usr/bin/env python3
import sys
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import axes3d, Axes3D
plt.rcParams['text.usetex'] = True
f = open( sys.argv[1], 'r' )
current_sigma = 0.0
current_rho = 0.0
current_beta = 0.0
sigma_lst = []
rho_lst = []
beta_lst = []
x_lst = []
y_lst = []
z_lst = []
x_data = []
y_data = []
z_data = []
parameters = []
data = {}
size = 0
for line in f:
line = line.strip()
a = line.split()
if not a:
continue
if a[ 0 ] == "#":
if x_lst:
parameters_tuple = ( current_sigma, current_rho, current_beta )
parameters.append( parameters_tuple )
data_tuple = ( np.array( x_lst ), np.array( y_lst ), np.array( z_lst ) )
data[ parameters_tuple ] = data_tuple
x_lst.clear()
y_lst.clear()
z_lst.clear()
current_sigma = float( a[ 2 ] )
current_rho = float( a[ 4 ] )
current_beta = float( a [ 6 ] )
if current_sigma not in sigma_lst:
sigma_lst.append( current_sigma )
if current_rho not in rho_lst:
rho_lst.append( current_rho )
if current_beta not in beta_lst:
beta_lst.append( current_beta )
else:
x_lst.append( float( a[ 0 ] ) )
y_lst.append( float( a[ 1 ] ) )
z_lst.append( float( a[ 2 ] ) )
parameters_tuple = ( current_sigma, current_rho, current_beta )
parameters.append( parameters_tuple )
data_tuple = ( np.array( x_lst ), np.array( y_lst ), np.array( z_lst ) )
data[ parameters_tuple ] = data_tuple
sigma_n = len( sigma_lst )
sigma_idx = 1
for sigma in sigma_lst:
rho_n = len( rho_lst )
beta_n = len( beta_lst )
fig, ax = plt.subplots( rho_n, beta_n, figsize=(8, 8), sharey=True, sharex=True )
fig.suptitle( f'$\sigma={sigma}$')
#ax = Axes3D(fig) dos not work with ax indexing
rho_idx = 0
beta_idx = 0
for rho in rho_lst:
for beta in beta_lst:
parameters_tuple = ( sigma, rho, beta )
data_tuple = data[ parameters_tuple ]
ax[ rho_idx, beta_idx ].plot( data_tuple[ 1 ], data_tuple[ 2 ], linewidth=1.0 )
if beta_idx == 0:
ax[ rho_idx, beta_idx ].set_ylabel( f'$\\rho={rho}$' )
if rho_idx == rho_n-1:
ax[ rho_idx, beta_idx ].set_xlabel( f'$\\beta={beta}$' )
beta_idx = beta_idx + 1
beta_idx = 0
rho_idx = rho_idx + 1
plt.savefig( f"{sys.argv[2]}-{sigma_idx}.png" )
sigma_idx = sigma_idx + 1
plt.close(fig)
......@@ -16,9 +16,9 @@ void solveParallelODEs( const char* file_name )
const Real output_time_step = 0.1;
const Real c_min = 1.0;
const Real c_max = 5.0;
const int c_vals = 11.0;
const int c_vals = 5.0;
const Real c_step = ( c_max - c_min ) / ( c_vals - 1 );
const int time_steps = final_t / tau + 1;
const int time_steps = final_t / output_time_step + 2;
Vector results( time_steps * c_vals, 0.0 );
auto results_view = results.getView();
......@@ -45,10 +45,12 @@ void solveParallelODEs( const char* file_name )
std::fstream file;
file.open( file_name, std::ios::out );
for( int k = 0; k < time_steps;k++ )
file << k * output_time_step << " ";
file << std::endl;
for( int i = 0; i < c_vals; i++ )
{
Real t = k * output_time_step;
file << t << " ";
for( int i = 0; i < c_vals; i++ )
for( int k = 0; k < time_steps;k++ )
file << results.getElement( k * c_vals + i ) << " ";
file << std::endl;
}
......@@ -56,8 +58,10 @@ void solveParallelODEs( const char* file_name )
int main( int argc, char* argv[] )
{
solveParallelODEs< TNL::Devices::Host >( "StaticODESolver-SineParallelExample-Host.out" );
TNL::String file_name( argv[ 1 ] );
file_name += "/StaticODESolver-SineParallelExample-result.out";
solveParallelODEs< TNL::Devices::Host >( file_name.getString() );
#ifdef HAVE_CUDA
solveParallelODEs< TNL::Devices::Cuda >( "StaticODESolver-SineParallelExample-Cuda.out" );
solveParallelODEs< TNL::Devices::Cuda >( file_name.getString() );
#endif
}
#!/usr/bin/env python3
import sys
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['text.usetex'] = True
f = open( sys.argv[1], 'r' )
data_lst = []
size = 0
for line in f:
line = line.strip()
a = line.split()
aux = []
for num in a:
aux.append( float( num ) )
data_lst.append( aux )
arrays = []
for a in data_lst:
arrays.append( np.array( a ) )
#x = arrays[ 0 ]
#arrays.remove( 0 )
n = len( arrays )
print( n )
fig, ax = plt.subplots( 1, n-1, figsize=(15, 3), sharey=True )
#fig, ax = plt.subplots( 1, n-1, sharey=True )
idx = 0
for array in arrays:
if idx > 0:
ax[ idx - 1 ].plot(arrays[0], array, linewidth=2.0)
ax[ idx - 1 ].set_xlabel( "t" )
ax[ idx - 1 ].set_ylabel( "u(t)" )
# ax[ idx - 1 ].axis('equal')
idx = idx + 1
plt.savefig( sys.argv[2] )
plt.close(fig)
......@@ -84,7 +84,7 @@ The first example solves ODE given by
\f[ u( 0 ) = 0, \f]
where \f$ c \f$ is a constant. We will solve it in parallel ODEs with different values \f$ c \in \langle c_{min}, c_{max} \rangle \f$. The code reads as follows:
where \f$ c \f$ is a constant. We will solve it in parallel ODEs with different values \f$ c \in \langle c_{min}, c_{max} \rangle \f$. The exact solution can be found [here](https://www.wolframalpha.com/input?i=y%27%28t%29+%3D++t+sin%28+a+t%29). The code reads as follows:
\includelineno Solvers/ODE/StaticODESolver-SineParallelExample.h
......@@ -94,7 +94,8 @@ Note, how we pass the value of parameter `c` to the lambda function `f`. The met
The result of this example looks as follows:
\include StaticODESolver-SineParallelExample-Host.out
\image{inline} html StaticODESolver-SineParallelExample.png ""
### Solving the Lorenz system in parallel
......@@ -125,6 +126,12 @@ It is very similar to the previous one. There are just the following changes:
The result looks as follows:
\image{inline} html StaticODESolver-LorenzParallelExample-1.png ""
\image{inline} html StaticODESolver-LorenzParallelExample-2.png ""
\image{inline} html StaticODESolver-LorenzParallelExample-3.png ""
\image{inline} html StaticODESolver-LorenzParallelExample-4.png ""
## Non-static ODE Solvers
In this section, we will show how to solve simple [heat equation](https://en.wikipedia.org/wiki/Heat_equation) in 1D. The heat equation is a parabolic partial differential equation which reads as follows:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment