//==============================================================================
//
//      Тестовая задача реализации параллельных вычислений OpenMP
//      (с) Притыкин Д. Е., к.т.н., доцент кафедры "Теоретическая механика"
//
//==============================================================================
#include        "openmp_train.h"
const   int     VEHICLES = 200;    
const   int     N = VEHICLES + 1;
const   double   c = 1.5e7;
double   dt = 1e-7;
double   t = 0;
double   tk = 10.0;
double  m[N]; // массы подвижных единиц (ПЕ)
double  y[2*N]; // фазовые координаты ПЕ
double  T[N+1]; // усилия в сцепках ПЕ
double  F[N]; // вектор внешних сил
double print_time = 0.0;
bool first_print = true;
timeval time0;
timeval time1;
double  calc_time = 0;
int i = 0;
int j = 0;
double dT = 0.01;
//------------------------------------------------------------------------------
//
//------------------------------------------------------------------------------
void OdeStep(double* mass, double* Y, double* Tc, double* force, double deltaT)
{
    #pragma omp for nowait
    for (i = 0; i < N; i++)
    {
        Tc[i+1] = c*(Y[i] - Y[i+1]);
        Tc[0] = 0.0;
        Tc[N] = 0.0;
        Y[i] = Y[i] + Y[i+N]*deltaT;
        Y[i+N] = Y[i+N] + deltaT*(force[i] + Tc[i] - Tc[i+1])/mass[i];
    }    
}
//------------------------------------------------------------------------------
//
//------------------------------------------------------------------------------
int main(void)
{
    // Инициализация
    m[0] = 129*1000;
    F[0] = 450*1000;
        
    for (i = 1; i < N; i++)
    {
        m[i] = 58*1000;
        F[i] = 0.0;
    }
    
    for (i = 0; i < 2*N; i++)
    {
        y[i] = 0.0;
    }
    
    for (i = 0; i < N+1; i++)
    {
        T[i] = 0.0;
    }   
      
    // Инициализация OpenMP
    omp_set_dynamic(0); // Запрещаем менять число потокв при выполнении
    omp_set_num_threads(8); // Задаем 8 потоков
          
    // Засекаем время начала расчета
    gettimeofday(&time0, NULL);  
// Задаем регион кода, вычисляемый параллелльно
#pragma omp parallel shared(m, y, T, F) private(i)
{
    while (t <= tk)
    {
        // Выполняем шаг интегрирования
        OdeStep(m, y, T, F, dt);  
        
        // Выводим результат на экран
       
        if ( ( t >= dT*j) || (j == 0) )
        {
                printf("%5i| time: %8.2f, с | progress: %5.1f % | T1: %8.1f, кН\n", j, t, t*100/tk, T[1]/1000.0);
                j++;            
        }
        
        
        t += dt;
    }
}
    
    // Засекаем время окончания расчета
    gettimeofday(&time1, NULL);
    
    // Вычисляем время расчета
    calc_time = ( time1.tv_sec + time1.tv_usec/1e6 ) - ( time0.tv_sec + time0.tv_usec/1e6 );
    
    printf(" time: %8.2f, с | progress: %5.1f % | T1: %8.1f, кН\n", t, t*100/tk, T[1]/1000.0);
    printf("\nCalculation time: %f, s\n", calc_time);    
    
    return 0;
}

#pragma         once
#include        <stdio.h>
#include        <math.h>
#include        <omp.h>
#include        <stdlib.h>
#include        <sys/time.h>