#include <iostream>
#include <fstream>
#include <math.h>
#include <omp.h>
using namespace std;
const int n = 100;
int i00 = n / 4;
int j00 = n / 4;
//определение источника
double diraca(int x)
{
if (x == 0)
return 1.0;
else
return 0.0;
}
double max(double a, double b)
{
if (a > b)
return a;
else
return b;
}
int main()
{
double V, U;
U = 1.0; V = 1.0;
double D = 0.01;
const double h = 100.0 / double(n);
double a1 = D / (h*h) + U / h;
double b1 = D / (h*h) + V / h;
double c1 = D / (h*h);
double d1 = D / (h*h);
double Q1 = 10.0;
int i, j;
double e1 = 2.0*D / (h*h) + U / h + 2.0*D / (h*h) + V / h;
double e2 = a1 + c1;
double e3 = b1 + d1;
double tau = 0.9;
double f1[n + 1][n + 1];
//решение на n-шаге
double C0[n + 1][n + 1];
//решение на (n+1/2)-шаге
double C1[n + 1][n + 1];
//решение на (n+1)-шаге
double C2[n + 1][n + 1];
//прогоночные коэфф.-ты
double P[n + 1], Q[n + 1];
//начальные условия распределения примеси
for (i = 0; i < n + 1; i++)
for (j = 0; j < n + 1; j++)
{
f1[i][j] = Q1 * diraca(i - i00) * diraca(j - j00);
C2[i][j] = 0.0;
C1[i][j] = 0.0;
C0[i][j] = 0.0;
}
double error = 1.0;
double t=omp_get_wtime();
while (error > 1e-7)
{
#pragma omp parallel for private(i,j) shared(C0,C2)
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
C0[i][j] = C2[i][j];
double a, b, c, d[n + 1];
//граничные условия шаг (n+1/2)
#pragma omp parallel for private(i) shared(C0,C1)
for (i = 0; i < n+1; i++)
{
C1[n][i] = C0[n - 1][i];
C1[i][n] = C0[i][n - 1];
C1[0][i] = 0.0;
C1[i][0] = 0.0;
}
//прогонка по столбцам
#pragma omp parallel for private(j) shared(a,b,c,i,d,C0,P,Q,C1)
for (j = 1; j < n; j++)
{
a = -a1; b = -(1 / tau + e2); c = -c1;
i = 1;
d[1] = d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
i = n-1;
d[n - 1] = b1 * C0[i][j-1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
for (i = 2; i < n-1; i++)
d[i] = b1 * C0[i][j - 1] + d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
P[1] = c / b;
Q[1] = -d[1] / b;
for (i = 2; i < n; i++)
{
P[i] = c / (-a * P[i - 1] + b);
Q[i] = (-d[i] + a * Q[i - 1]) / (-a * P[i - 1] + b);
}
C1[n - 1][j] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
for (i = n - 2; i >= 1; i--)
{
C1[i][j] = P[i] * C1[i + 1][j] + Q[i];
}
}
//граничные условия шаг (n+2)
#pragma omp parallel for private(i) shared(C1,C2)
for (i = 0; i < n+1; i++)
{ C2[n][i] = C1[n - 1][i];
C2[i][n] = C1[i][n - 1];
C2[0][i] = 0.0;
C2[i][0] = 0.0;
}
//прогонка по строкам
#pragma omp parallel for private(i) shared(a,b,c,j,d,C1,P,Q,C2)
for (i = 1; i < n; i++)
{
a = -b1; b = -(1 / tau + e3); c = -d1;
j = 1;
d[1] = C1[i+1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
j = n-1;
d[n - 1] = a1 * C1[i - 1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
for (j = 2; j < n-1; j++)
d[j] = a1 * C1[i - 1][j] + c1 * C1[i + 1][j] - (1 / tau - e2)*C1[i][j] + f1[i][j];
P[1] = c / b;
Q[1] = -d[1] / b;
for (j = 2; j < n; j++)
{
P[j] = c / (-a * P[j - 1] + b);
Q[j] = (-d[j] + a * Q[j - 1]) / (-a * P[j - 1] + b);
}
C2[i][n - 1] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
for (j = n - 2; j >= 1; j--)
{
C2[i][j] = P[j] * C2[i][j + 1] + Q[j];
}
}
error = 0.0;
#pragma omp parallel for private(i,j) reduction(+:error)
for (i = 1; i < n; i++)
for (j = 1; j < n; j++)
error += pow(fabs(C2[i][j] - C0[i][j]), 2);
error = pow(error / (n*n), 0.5);
}
t=omp_get_wtime()-t;
FILE *f;
f=fopen("result.txt","w");
for(i=1;i<n;i++)
for(j=1;j<n;j++)
fprintf(f,"%i %i %lf\n",i,j,C2[i][j]);
fclose(f);
printf("Time=%lf\n",t);
return 0;
}
