---

#define _POSIX_C_SOURCE 199309L

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define NX 1000
#define NY 1000
#define STEPS 500
#define L 0.01
#define DX (L/NX)
#define DY (L/NY)
#define Kk 6.196522218e-6
#define DT (0.9 * DX * DY / (4 * Kk))

#define PI 3.141592653589793

// Вспомогательные функции
__host__ __device__ double lamda(double u) {
    return 1.45 + 0.023 * u - 2e-6 * u * u;
}

__host__ __device__ double c(double u) {
    return 1 / (2.25 * 0.001 - 6.08e-10 * u * u);
}

__host__ __device__ double p(double u) {
    return 7860.0 + 41500.0 / u;
}

__host__ __device__ double u0(double x, double y) {
    x /= L;
    y /= L;
    return 400.0 + 50.0 * x * x + 50.0 * y * y;
}

__host__ __device__ double m1(double y, double t) {
    y /= L;
    t /= (DT * STEPS);
    return 400.0 + 50.0 * y * y + 100.0 * t * sqrt(16.0 + 20.0 * y * y);
}

__host__ __device__ double m2(double y, double t) {
    y /= L;
    t /= (DT * STEPS);
    return 450.0 + 50.0 * y * y + 150.0 * t * (1.0 - y);
}

__host__ __device__ double m3(double x, double t) {
    x /= L;
    t /= (DT * STEPS);
    return 400.0 + 50.0 * x * x + 400.0 * t * cos(2.0 * PI * x) - 250.0 * t * x;
}

__host__ __device__ double m4(double x, double t) {
    x /= L;
    t /= (DT * STEPS);
    return 450.0 + 50.0 * x * x + 600.0 * t * cos(PI * x / 2.0);
}

// Ядро CUDA для обновления внутренних точек
__global__ void update_kernel(double* u_curr, double* u_next, int nx, int ny, double dx, double dy, double dt) {
    int i = blockIdx.x * blockDim.x + threadIdx.x + 1;
    int j = blockIdx.y * blockDim.y + threadIdx.y + 1;
    
    if (i < nx - 1 && j < ny - 1) {
        int idx = i * ny + j;
        
        double u_center = u_curr[idx];
        
        double u_left = u_curr[(i-1) * ny + j];
        double u_right = u_curr[(i+1) * ny + j];
        double u_bottom = u_curr[i * ny + (j-1)];
        double u_top = u_curr[i * ny + (j+1)];
        
        double lamda_left = lamda((u_left + u_center) / 2.0);
        double lamda_right = lamda((u_right + u_center) / 2.0);
        double lamda_bottom = lamda((u_bottom + u_center) / 2.0);
        double lamda_top = lamda((u_top + u_center) / 2.0);
        
        double flux_x = lamda_right * (u_right - u_center) - lamda_left * (u_center - u_left);
        double flux_y = lamda_top * (u_top - u_center) - lamda_bottom * (u_center - u_bottom);
        
        double div = c(u_center) * p(u_center);
        
        u_next[idx] = u_center + dt * (flux_x / (dx * dx) + flux_y / (dy * dy)) / div;
    }
}

// Ядро для установки граничных условий
__global__ void boundary_kernel(double* u, int nx, int ny, double dx, double dy, double t) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    
    if (i == 0 && j < ny) {
        u[0 * ny + j] = m1(j * dy, t);
    }
    else if (i == nx - 1 && j < ny) {
        u[(nx - 1) * ny + j] = m2(j * dy, t);
    }
    else if (j == 0 && i < nx) {
        u[i * ny + 0] = m3(i * dx, t);
    }
    else if (j == ny - 1 && i < nx) {
        u[i * ny + (ny - 1)] = m4(i * dx, t);
    }
}

// Функция замера времени
double get_time() {
    struct timespec ts;
    clock_gettime(CLOCK_MONOTONIC, &ts);
    return ts.tv_sec + ts.tv_nsec * 1e-9;
}

int main() {
    // Выделение памяти на хосте
    double *h_u_curr = (double*)malloc(NX * NY * sizeof(double));
    double *h_u_next = (double*)malloc(NX * NY * sizeof(double));
    
    if (!h_u_curr || !h_u_next) {
        printf("Ошибка выделения памяти на хосте!\n");
        return 1;
    }
    
    // Инициализация начальными условиями
    for (int i = 0; i < NX; i++) {
        for (int j = 0; j < NY; j++) {
            int idx = i * NY + j;
            h_u_curr[idx] = u0(i * DX, j * DY);
            h_u_next[idx] = 0.0;
        }
    }
    
    // Выделение памяти на устройстве
    double *d_u_curr, *d_u_next;
    cudaMalloc(&d_u_curr, NX * NY * sizeof(double));
    cudaMalloc(&d_u_next, NX * NY * sizeof(double));
    
    // Копирование начальных данных на GPU
    cudaMemcpy(d_u_curr, h_u_curr, NX * NY * sizeof(double), cudaMemcpyHostToDevice);
    cudaMemcpy(d_u_next, h_u_next, NX * NY * sizeof(double), cudaMemcpyHostToDevice);
    
    // Настройка сетки блоков и потоков
    dim3 threads(16, 16);
    dim3 blocks((NX - 2 + threads.x - 1) / threads.x, (NY - 2 + threads.y - 1) / threads.y);
    dim3 blocks_boundary((NX + threads.x - 1) / threads.x, (NY + threads.y - 1) / threads.y);
    
    double start_time = get_time();
    
    // Основной цикл по времени
    for (int step = 0; step < STEPS; step++) {
        double t = DT * (step + 1);
        
        update_kernel<<<blocks, threads>>>(d_u_curr, d_u_next, NX, NY, DX, DY, DT);
        
        boundary_kernel<<<blocks_boundary, threads>>>(d_u_next, NX, NY, DX, DY, t);
        
        double *temp = d_u_curr;
        d_u_curr = d_u_next;
        d_u_next = temp;
    }
    
    cudaDeviceSynchronize();
    double end_time = get_time();
    
    // Копирование результата обратно на хост
    cudaMemcpy(h_u_curr, d_u_curr, NX * NY * sizeof(double), cudaMemcpyDeviceToHost);
    
    printf("\nРезультат:\n");
    printf("Температура в центре: %.2f\n", h_u_curr[(NX/2) * NY + NY/2]);
    printf("Время расчёта: %.3f секунд\n", end_time - start_time);
    
    // Освобождение памяти
    cudaFree(d_u_curr);
    cudaFree(d_u_next);
    free(h_u_curr);
    free(h_u_next);
    
    return 0;
}