This file provides a brief documentation and information related to the fourth (last) Homework of the course "Scientific Programming for Engineers", fall 2019.
This homework is done by O. Ashtari and A. Sieber.
Last update: 15.01.2020
The first part of the Homework aims at using the Pybind11 library in order to create Python binding for the C++ particles code. The ultimate goal of this part consists in allowing the user to run this code through a Python interface.
The second part of the Homework focuses on the planet trajectory simulation module of the particles code. There, a Python optimization routine is written which main purpose is to modify the initial planets states (one planet at a time) in order to minimize the error of their computed trajectories when compared to reference trajectories.
After cloning this repository on your computer, you should build the program in a desired directory based on the CMakeLists.txt file in the root folder. To do so, from within the destination directory, you can build the executable file(s) using command line below:
$ cmake -DCMAKE_PREFIX_PATH=<include> -DCMAKE_LIBRARY_PATH=<library> -DUSE_FFT=<status> -DUSE_PYTHON=<status> <CMakeLists> $ make -j
where the command cmake is basically followed by the address of the directory where the CMakeLists.txt is located (denoted by <CMakeLists>.) Additional options are available as follows:
- -DCMAKE_PREFIX_PATH specifies the path where you want the program to look for headers (.h files). <include> in the command above will be replaced by the address.
- -DCMAKE_LIBRARY_PATH specifies the path where you want the program to look for libraries (.a files). <library> in the command above will be replaced by the address.
- -DUSE_FFT determines whether the program should look for FFTW library/header in the specified directories or not. <status> can either be ON (default) or OFF.
- -DUSE_PYTHON determines whether the program should look for pybind11 library/header in the specified directories or not. <status> can either be ON (default) or OFF.
1- The minimum requirement for CMake is the version 3.1.
2- googletest must be either installed in a default/specified directory or cloned in the lib folder to be build the executables. The library can be cloned from https://github.com/google/googletest.
3- pybind11 must be installed in a default directory where the program looks for headers/libraries, otherwise you should provide the directory when using cmake command (see explanation for the sample cmake command line.) The library can be cloned from https://github.com/pybind/pybind11.
Executable files will be built in <build_directory>/src. In this program we only focus on using python scripts main.py and error_minimization.py.
The optimization routine error_minimization.py is intended to minimize the error on Mercury trajectory (or any other planet of the solar system) by scaling the initial velocity of the planet of interest. The minimization is performed by computing the error between the computed planet trajectory and a reference planet trajectory. This routine calls main.py, therefore, in order to run it, the user must first make sure she/he can launch the particles code through the python interface. Once the requirement is satisfied the error_minimization.py routine can be called as follow:
$ python3 error_minimization.py directory_comp directory_ref planet input_file scale nb_steps freq
- directory_comp, the path to the directory containing the computed trajectory data.
- directory_ref, the path to the directory containing the reference trajectory data.
- planet, the planet of interest.
- input_file, the file containing the initial state of the different planets (the path to the file must be included).
- scale, the initial scaling factor (initial guess of the minimization).
- nb_steps, the number of time steps of the simulation within the particles code.
- freq, the dumping frequency at which the particles code writes outputs.
As an example a minimization on Mercury trajectory over 365 days with a dumping frequency of 1 day could be run as such:
$ python3 error_minimization.py dumps trajectories mercury init.csv 1.0 365 1
The error_minimization.py routine then prints the scaling factor minimizing the error, the value of this error and the amount of minimization iterations needed to reach an optimum. It moreover plots the evolution of the error versus the scaling factor as shown in the link below:
The above graph shows that the optimized scaling is evaluated at 0.39895. Comparing this result to the exact scaling of 0.4 (found by dividing the reference initial velocity by the one stored in init.csv) leads to a relative error of less than 0.26% in the optimization procedure.
To comment on this function overload, let's first discuss what is the role of createComputes in the whole code.
In the ParticlesFactoryInterface class, createComputes is a Real-to-voind function which is to be defined later, if needed. In the constructor of MaterialPointsFactory, createComputes is defined to be the default function (i.e. to be createDefaultComputes.) In the default function, ComputeTemperature is added to the system evolution object. Similarly in the PlanetsFactory where createComputes is defined to be createDefaultComputes in which verlet is added to the system evolution object. As a summary, by constructing an object of type MaterialPointsFactory or PlanetsFactory, createComputes takes a default definition. It is then used inside the first definition of createSimulation method whose arguments are file name and time step size.
createSimulation is overloaded such that createComputes is redefined by receiving and substituting the python functor with the desired functionality. The python object (functor) includes initializations and adding compute objects explained above. The whole idea behind overloading createSimulation and using the templated functor is to give the flexibility to the python user to easily tailor computing cobjects (i.e. ComputeTemperature or verlet etc.) based on his/her needs. The reason is that when a C++ code is wrapped to be used through python, the C++ side is supposed to not be touched anymore. However, default settings for material point, for instance, is hard-coded inside its factory. Therefore, to change default values, the python user should manipulate the C++ files. Or assume, like the previous homework, one wants to get length, conductivity, etc. from command line arguments. Or one wants to develop the code such that these variables are read from a file. Using the overloaded createSimulation with having very little insight into how the C++ code is designed, one can do all these by just writing few lines in python, with neither changing any C++ file nor editing the binding file.