From 2de388af0dbbe410c9c7d2897f8cc8b25da9ad5a Mon Sep 17 00:00:00 2001
From: eug-vs <eugene@eug-vs.xyz>
Date: Mon, 21 Mar 2022 23:04:40 +0300
Subject: feat: add source code

---
 src/report/report.tex |  12 ++
 src/source.cpp        | 501 ++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 513 insertions(+)
 create mode 100755 src/source.cpp

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;
+}
+
+
+
-- 
cgit v1.2.3