feat: add source code

2 files changed, 513 insertions, 0 deletions

diff --git a/src/report/report.tex b/src/report/report.tex
index e4cd7bb..aef5cfa 100644
--- a/src/report/report.tex
+++ b/src/report/report.tex
@@ -229,4 +229,16 @@
     \caption{Вихревой след за ступенькой}
 \end{figure}
 
+\chapter{Код программы}
+\lstinputlisting[language=C, caption=Исходный код, label=lst:code]{./src/source.cpp}
+
+\chapter{Список использованных источников}
+\begin{enumerate}
+  \item{Patankar SV. Numerical Heat Transfer and Fluid Flow, 1984}
+  \item{Versteeg HK, Malalasekera W. An Introduction to Computational Fluid Dynamics, 1995}
+  \item{Devic I. Fluid simulation with SIMPLE method using graphic processors, 2014}
+  \item{Ferziger JH, Peric M. Computational Methods for Fluid Dynamics, 2002}
+  \item{Андерсон Д, Таннехилл Дж, Плетчер Р. Вычислительная гидромеханика и теплообмен, 1990}
+\end{enumerate}
+
 \end{document}
diff --git a/src/source.cpp b/src/source.cpp
new file mode 100755
index 0000000..3c35588
--- /dev/null
+++ b/src/source.cpp
@@ -0,0 +1,501 @@
+// This is not my code.

+// I have no idea how it works.

+// It probably requires paid software.

+// I'm planning my own implementation in Python.

+#include <stdlib.h>

+#include <stdio.h>

+#include <time.h>

+#include <cmath>

+#include "TECIO.h"

+

+#ifndef NULL

+#define NULL 0

+#endif

+

+enum FileType { FULL = 0, GRID = 1, SOLUTION = 2 };

+

+#define BD 20

+#define BX 50

+#define BY 10

+#define NX 2000

+#define NY 200

+#define KT 10

+#define AI -2./3.

+#define DT 0.00001

+#define RX 100

+#define RY 40

+

+#define IY 40

+#define IX 80

+

+#define Re 400

+#define GR 0.0001

+#define T 6000

+

+double dxx[NX], dyy[NY], xi[NX][NY], yi[NX][NY], a[NX*NY], b[NX*NY], c[NX*NY];

+

+double dx[NX][NY], dy[NX][NY];

+

+

+

+

+void PRINT_PLT(void) {

+	double SolTime;

+	INTEGER4 Debug, I, J, III, DIsDouble, VIsDouble, IMax, JMax, KMax, ZoneType, StrandID, ParentZn, IsBlock;

+	INTEGER4 ICellMax, JCellMax, KCellMax, NFConns, FNMode, ShrConn, FileType;

+

+	Debug = 1;

+	VIsDouble = 1;

+	DIsDouble = 1;

+	IMax = NY;

+	JMax = NX;

+	KMax = 1;

+	ZoneType = 0;      /* Ordered */

+	SolTime = 360.0;

+	StrandID = 0;     /* StaticZone */

+	ParentZn = 0;      /* No Parent */

+	IsBlock = 1;      /* Block */

+	ICellMax = 0;

+	JCellMax = 0;

+	KCellMax = 0;

+	NFConns = 0;

+	FNMode = 0;

+	ShrConn = 0;

+	FileType = FULL;

+

+	/*

+	* Open the file and write the tecplot datafile

+	* header information

+	*/

+	I = TECINI112((char*)"SIMPLE DATASET",

+		(char*)"X Y U V",

+		(char*)"UV.plt",

+		(char*)".",

+		&FileType,

+		&Debug,

+		&VIsDouble);

+

+	/*

+	* Write the zone header information.

+	*/

+	I = TECZNE112((char*)"Simple Zone",

+		&ZoneType,

+		&IMax,

+		&JMax,

+		&KMax,

+		&ICellMax,

+		&JCellMax,

+		&KCellMax,

+		&SolTime,

+		&StrandID,

+		&ParentZn,

+		&IsBlock,

+		&NFConns,

+		&FNMode,

+		0,              /* TotalNumFaceNodes */

+		0,              /* NumConnectedBoundaryFaces */

+		0,              /* TotalNumBoundaryConnections */

+		NULL,           /* PassiveVarList */

+		NULL,           /* ValueLocation = Nodal */

+		NULL,           /* SharVarFromZone */

+		&ShrConn);

+

+	/*

+	* Write out the field data.

+	*/

+	III = IMax * JMax;

+	I = TECDAT112(&III, &dx[0][0], &DIsDouble);

+	I = TECDAT112(&III, &dy[0][0], &DIsDouble);

+	I = TECDAT112(&III, &xi[0][0], &DIsDouble);

+	I = TECDAT112(&III, &yi[0][0], &DIsDouble);

+

+	I = TECEND112();

+

+

+}

+

+

+void PRINT_VECT_VTK(double U[NX][NY], double V[NX][NY], double *x, double *y) {

+	int i, j;

+

+	FILE *vect_field;

+

+	/* File grid.vtk*/

+	if ((vect_field = fopen("a.vtk", "w")) == NULL)

+		printf("The grid file could not be open\n");

+	else{

+		fprintf(vect_field, "%d\n", NX);

+		fprintf(vect_field, "%d\n", NY);

+

+		for (i = 0; i <= NX - 1; i++)

+			fprintf(vect_field, "%f\n", x[i]);

+

+		for (j = 0; j <= NY - 1; j++)

+			fprintf(vect_field, "%f\n", y[j]);

+

+		for (i = 0; i <= NY - 1; i++)

+			for (j = 0; j <= NX - 1; j++){

+			fprintf(vect_field, "%f\t %f\t\n", U[j][i], V[j][i]);

+			}

+	}

+	fclose(vect_field);

+

+}

+

+void PRINT_VTK(double U[NX][NY], double V[NX][NY], double *x, double *y) {

+	int i, j;

+

+	FILE *vect_field;

+

+

+	/* File grid.vtk*/

+	if ((vect_field = fopen("b.vtk", "w")) == NULL)

+		printf("The grid file could not be open\n");

+	else{

+		fprintf(vect_field, "# vtk DataFile Version 2.0\n");

+		fprintf(vect_field, "Sample rectilinear grid\n");

+		fprintf(vect_field, "ASCII\n");

+		fprintf(vect_field, "DATASET RECTILINEAR_GRID\n");

+		fprintf(vect_field, "DIMENSIONS %d %d \n", NX, NY);

+		fprintf(vect_field, "X_COORDINATES %d double\n", NX);

+		for (i = 0; i <= NX - 1; i++)

+			fprintf(vect_field, "%f\n", x[i]);

+		fprintf(vect_field, "Y_COORDINATES %d double\n", NY);

+		for (j = 0; j <= NY - 1; j++)

+			fprintf(vect_field, "%f\n", y[j]);

+//		fprintf(vect_field, "Z_COORDINATES %d double\n", nz);

+//		fprintf(vect_field, "%f\n", z);

+		fprintf(vect_field, "POINT_DATA %d\n", (NX)*(NY));

+		fprintf(vect_field, "VECTORS vectors double\n");

+		for (i = 0; i <= NY - 1; i++)

+			for (j = 0; j <= NX - 1; j++){

+				fprintf(vect_field, "%f\t %f\t %f\n", U[j][i], V[j][i], 0.0);

+			}

+	}

+	fclose(vect_field);

+

+}

+

+

+void inixh(double *a)

+{

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+

+			k = i *NY + j;

+			if (i<RX && j>RY)

+				a[k] = 1.;

+			else

+				a[k] = 0.;

+		}

+	}

+}

+void inih(double *a)

+{

+	int i, j, k;

+

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+			a[k] = 0.;

+		}

+	}

+}

+void USh(double *a, double *b, double *c, double *d)

+{

+	double u1;

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+

+			if (i < NX - 2 && i>0 && j < NY - 1 && j>0 && ((i <= RX && j <= RY) != 1)){

+				u1 = -((pow(a[k + NY], 2) - pow(a[k - NY], 2))*NX / 2.0f + (a[k + 1] * 0.5*(b[k] + b[k + NY]) - a[k - 1] * 0.5*(b[k - 1] + b[k + NY - 1]))*NY*KT / 2.0);

+				u1 = u1 + (1.0 / Re)*((a[k + NY] - 2.0*a[k] + a[k - NY])*NX*NX + (a[k + 1] - 2.0*a[k] + a[k - 1])*NY*NY*KT*KT);

+				u1 = a[k] + DT*(u1 - 1.0*NX*(c[k + NY] - c[k]));

+				d[k] = u1;

+			}

+			else

+				d[k] = a[k];

+

+		}

+	}

+}

+void VSh(double *a, double *b, double *c, double *d)

+{

+

+	double v1;

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+

+			if (i < NX - 1 && i>0 && j < NY - 2 && j>0 && ((i <= RX && j <= RY) != 1)){

+				v1 = -((b[k + NY] * 0.5*(a[k] + a[k + 1]) - b[k - NY] * 0.5*(a[k - NY] + a[k - NY + 1]))*NX / 2.0 + (b[k + 1] * b[k + 1] - b[k - 1] * b[k - 1])*KT*NY / 2.0);

+				v1 = v1 + (1.0f / Re)*((b[k + NY] - 2.0f*b[k] + b[k - NY])*NX*NX + (b[k + 1] - 2.0f*b[k] + b[k - 1])*NY*NY*KT*KT);

+				v1 = b[k] + DT*(v1 - (1.0f*NY*KT)*(c[k + 1] - c[k]));

+				d[k] = v1;

+			}

+			else

+				d[k] = b[k];

+		}

+	}

+}

+void sumah(double *a, double *b)

+{

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+			b[k] += a[k];

+		}

+	}

+}

+void rubnixh(double *a)

+{

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+

+			if (i == 0 && j>RY)

+				a[k] = 1.;

+			if (i == NX - 2)

+				a[k] = a[k - NY];

+			if (j == 0)

+				a[k] = 0.;

+			if (j == NY - 1)

+				a[k] = a[k - 1];

+			if (j == RY && i < RX)

+				a[k] = 0.;

+			if (i == RX && j <= RY)

+				a[k] = a[k + NY] * AI;

+		}

+	}

+}

+void rubniyh(double *a)

+{

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+

+			if (j == 0)

+				a[k] = AI*a[k + 1];

+			if (j == NY - 2)

+				a[k] = a[k - 1];

+			if (i == 0)

+				a[k] = 0.;

+			if (i == NX - 1)

+				a[k] = a[k - NY];

+			if (j == RY && i <= RX)

+				a[k] = AI*a[k + 1];

+			if (i == RX && j < RY)

+				a[k] = 0.;

+		}

+	}

+}

+void rubniph(double *a)

+{

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+

+

+			if (j == 0 && i != 0 && i != NX - 1)

+				a[k] = a[k + 1];

+			if (j == NY - 1 && i != 0 && i != NX - 1)

+				a[k] = 0.0;

+			if (i == 0 && j != 0 && j != NY - 1)

+				a[k] = 0;

+			if (i == NX - 1 && j != 0 && j != NY - 1)

+				a[k] = 0.;

+			if (i == RX && j <= RY)

+				a[k] = a[k + NY];

+			if (i <= RX && j == RY)

+				a[k] = a[k + 1];

+		}

+	}

+}

+void korekh(double *a, double *b, double *c)

+{

+	int i, j, k;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+

+			if (i>0 && i < NX - 2 && ((i <= RX && j <= RY) != 1)){

+				a[k] -= DT*NX*(c[k + NY] - c[k]);

+			}

+

+			if (j>0 && j < NY - 2 && ((i <= RX && j <= RY) != 1)){

+				b[k] -= DT*NY*KT*(c[k + 1] - c[k]);

+			}

+

+		}

+	}

+}

+

+void poissh(double *a, double *b, double *c, double *d)

+{

+	int i, j, k;

+	double corr;

+	for (i = 0; i < NX; i++) {

+		for (j = 0; j < NY; j++) {

+			k = i *NY + j;

+			corr = c[k];

+			if (i > 0 && j > 0 && i < NX - 1 && j < NY - 1 && ((i <= RX && j <= RY) != 1)){

+				c[k] = -(1.0 / (2 * (DT*NX*NX + DT*NY*NY*KT*KT)))*((-DT*NX*NX)*c[k + NY] + (-DT*NX*NX)*c[k - NY] + (-DT*NY*NY*KT*KT)*c[k + 1] + (-DT*NY*NY*KT*KT)*c[k - 1] + (1.0*NX)*(a[k] - a[k - NY]) + (1.0*NY*KT)*(b[k] - b[k - 1]));

+			}

+			corr -= c[k];

+			if (corr < 0)

+				corr *= (-1);

+			d[k] = corr;

+		}

+	}

+}

+

+

+

+int main()

+{

+

+	clock_t t2, t1;

+	FILE *fp, *fp1, *fp2;

+

+	double *a_h, *a_h1, *b_h, *b_h1, *c_h, *c_h1, *corr_h;

+

+	double  dx1, dy1;

+	int j, i, t, k;

+

+	fp = fopen("x.dat", "w");

+	fp1 = fopen("y.dat", "w");

+	fp2 = fopen("tlak.dat", "w");

+	t1 = clock();

+

+	a_h = (double*)malloc(sizeof(double)*NX*NY);

+	a_h1 = (double*)malloc(sizeof(double)*NX*NY);

+	b_h = (double*)malloc(sizeof(double)*NX*NY);

+	b_h1 = (double*)malloc(sizeof(double)*NX*NY);

+	c_h = (double*)malloc(sizeof(double)*NX*NY);

+	c_h1 = (double*)malloc(sizeof(double)*NX*NY);

+	corr_h = (double*)malloc(sizeof(double)*NX*NY);

+

+	inixh(a_h);

+	inixh(a_h1);

+	inih(b_h);

+	inih(b_h1);

+	inih(c_h);

+	inih(c_h1);

+

+

+	printf("Start");

+	for (t = 0; t<T; t++)

+	{

+

+		USh(a_h, b_h, c_h, a_h1);

+		VSh(a_h, b_h, c_h, b_h1);

+		rubnixh(a_h1);

+		rubniyh(b_h1);

+

+		for (i = 0; i<30; i++)

+			poissh(a_h1, b_h1, c_h1, corr_h);

+

+		rubniph(c_h1);

+		korekh(a_h1, b_h1, c_h1);

+		rubnixh(a_h1);

+		rubniyh(b_h1);

+		sumah(c_h1, c_h);

+		inih(c_h1);

+		rubniph(c_h);

+		USh(a_h1, b_h, c_h, a_h);

+		VSh(a_h, b_h1, c_h, b_h);

+		rubnixh(a_h);

+		rubniyh(b_h);

+

+		for (i = 0; i<30; i++)

+	    poissh(a_h, b_h, c_h1, corr_h);

+

+

+		rubniph(c_h1);

+		korekh(a_h, b_h, c_h1);

+		rubnixh(a_h);

+		rubniyh(b_h);

+		sumah(c_h1, c_h);

+		inih(c_h1);

+		rubniph(c_h);

+

+

+	}

+	for (i = 0; i<NX; i++)

+		for (j = 0; j<NY; j++) if (i != NX - 1)

+			fprintf(fp, "\n%d %d %f", i, j, a_h[i*NY + j]);

+	for (i = 0; i<NX; i++)

+		for (j = 0; j<NY; j++) if (j != NY - 1)

+			fprintf(fp1, "\n%d %d %f", i, j, b_h[i*NY + j]);

+	for (i = 0; i<NX; i++)

+		for (j = 0; j<NY; j++)

+			fprintf(fp2, "\n%d %d %f", i, j, c_h[i*NY + j]);

+	t2 = clock();

+	printf(" %f\n", ((double)(t2 - t1)) / CLOCKS_PER_SEC);

+

+	dx1 = 1.0*NX - 1;

+	dy1 = KT*NY - 1;

+	for (i = 0; i<NX; i++)

+		dxx[i] = 1.0*i / dx1;

+	for (i = 0; i<NY; i++)

+		dyy[i] = 1.0*i / dy1;

+	for (i = 0; i<NX; i++)

+	{

+		xi[i][0] = 0.0;

+		yi[i][0] = 0.0;

+		xi[i][NY - 1] = 0.0;

+		yi[i][NY - 1] = 0.0;

+	}

+	for (i = 0; i<NY; i++)

+	{

+		xi[0][i] = 1.0;

+		yi[0][i] = 0.0;

+		xi[NX - 1][i] = 0.0;

+		yi[NX - 1][i] = 0.0;

+	}

+	for (i = 1; i <= NX - 2; i++)

+		for (j = 1; j <= NY - 1; j++)

+		{

+		k = (i - 1)*NY + j;

+

+		xi[i][j] = 0.5*(a_h[k] + a_h[k + NY]);

+		}

+	for (i = 1; i <= NX - 1; i++)

+		for (j = 1; j <= NY - 2; j++)

+		{

+		k = i*NY + j - 1;

+

+		yi[i][j] = 0.5*(b_h[k] + b_h[k + 1]);

+		}

+

+	PRINT_VECT_VTK(xi, yi, dxx, dyy);

+	PRINT_VTK(xi, yi, dxx, dyy);

+

+	for (j = 0; j < NX; j++)

+	for (i = 0; i < NY ; i++)

+		 {

+			dx[j][i] = (1.0*j / dx1);

+			dy[j][i]=  (1.0*i / dy1);

+		}

+

+

+

+	PRINT_PLT();

+

+

+	printf("END");

+

+

+    return 0;

+}

+

+

+