1

code/CMakeLists.txt

@@ -0,0 +1,21 @@
+ PROJECT(main)
+ CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
+include(CheckCXXCompilerFlag)
+CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11)
+CHECK_CXX_COMPILER_FLAG("-std=c++0x" COMPILER_SUPPORTS_CXX0X)
+if(COMPILER_SUPPORTS_CXX11)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+elseif(COMPILER_SUPPORTS_CXX0X)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x")
+else()
+        message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11
+support. Please use a different C++ compiler.")
+endif() 
+SET(CMAKE_CXX_FLAGS_RELEASE "$ENV{CXXFLAGS} -O3 -Wall")
+find_package(OpenMP)
+if (OPENMP_FOUND)
+    set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+    set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif()
+AUX_SOURCE_DIRECTORY(. DIR_SRCS)
+ ADD_EXECUTABLE(main ${DIR_SRCS})

code/accuracy_test.cpp

@@ -0,0 +1,524 @@
+#include"accuracy_test.h"
+
+void accuracy_sinwave_1d()
+{
+	int mesh_set = 3; //用一系列网格尺寸等比例减小的网格来测试格式的精度
+	int mesh_number_start = 10; //最粗网格数
+	double length = 2.0; //计算域总长度为2，[0,2]
+	double CFL = 0.5; //CFL数
+
+	// 注意，为了防止时间推进和空间离散的精度不一致，dt = CFL * （dx^dt_ratio），这一点的原因也可以参看
+	// X. Ji et al. Journal of Computational Physics 356 (2018) 150–173 的 subsection 4.1
+	double dt_ratio = 1.0;
+	//end
+
+	double mesh_size_start = length / mesh_number_start; //均匀网格大小
+
+	//各套网格数量-尺寸-误差 开辟数组存一下
+	int* mesh_number = new int[mesh_set];
+	double* mesh_size = new double[mesh_set];
+	double** error = new double* [mesh_set];
+	//end
+
+	for (int i = 0; i < mesh_set; i++)
+	{
+		error[i] = new double[3]; //每一套网格都计算 L1 L2 L_inf 误差
+		mesh_size[i] = mesh_size_start / pow(2, i); //等比变化
+		mesh_number[i] = mesh_number_start * pow(2, i); //等比变化
+
+		//计算每一套网格的工况,得到对应的误差
+		accuracy_sinwave_1d(CFL, dt_ratio, mesh_number[i], error[i]);
+		//end
+	}
+
+	output_error_form(CFL, dt_ratio, mesh_set, mesh_number, error);
+}
+void accuracy_sinwave_1d(double& CFL, double& dt_ratio, int& mesh_number, double* error)
+{
+	Runtime runtime;//记录一下计算时间
+	runtime.start_initial = clock();
+
+	Block1d block; //存储一些一维算例相关的基本信息
+
+	block.nodex = mesh_number;
+	block.ghost = 3; // ghost cell 数。一般来讲 二阶-三阶 需要两个，四阶-五阶 需要三个。具体算法则需具体讨论
+
+	double tstop = 2.0; //计算一个sinwave的周期
+
+	block.CFL = CFL; // CFL数
+
+
+
+	K = 4; //内部自由度
+	Gamma = 1.4; //比热容是由内部自由度决定的，这里为了简单，直接手动给出
+
+	tau_type = Euler; //无粘还是粘性问题
+
+	//准备通量函数
+	//函数指针 flux_function 可以方便的切换 使用的通量求解器
+	//如 L-F， HLLC, GKS
+	flux_function = GKS; //GKS
+	//通常 GKS是指 2阶gks，但通过一些近似，可以变形为
+	//一阶GKS（gks1st） 一阶KFVS(kfvs1st) 等，在gks1dsolver中枚举出来
+	gks1dsolver = gks2nd; //现在选择二阶gks
+	c1_euler = 0.0; //仅对于无粘流动，gks里面的人工碰撞时间系数c1
+	c2_euler = 0.0; //对于无粘或者粘性流动，gks里面的人工碰撞时间系数c2
+	//end
+
+	//边界条件函数设定
+	Fluid1d* bcvalue = new Fluid1d[2]; //一维左右dirichlet边界的值
+	BoundaryCondition leftboundary(0);
+	BoundaryCondition rightboundary(0);
+	leftboundary = periodic_boundary_left;
+	rightboundary = periodic_boundary_right;
+	//end
+
+	//设置重构函数
+	cellreconstruction = WENO5_AO; //以单位为中心的重构，重构出界面左值和右值
+	wenotype = wenoz; //非线性权的种类,linear就意味着采用线性权
+	reconstruction_variable = characteristic; //重构变量的类型，分为特征量和守恒量
+	g0reconstruction = Center_collision;//平衡态函数g0的重构方式（传统黎曼求解器用此函数得到粘性项所需要的一阶导数）
+	//end
+
+	//显式时间推进选取，支持N步M阶时间精度，N目前写死了小于等于5
+	timecoe_list = S2O4;//这叫做2步4阶 two-stage fourth-order
+	Initial_stages(block); //赋值一下
+	//end
+
+	// allocate memory for 1-D fluid field
+	// in a standard finite element method, we have 
+	// first the cell average value, N
+	block.nx = block.nodex + 2 * block.ghost;
+	block.nodex_begin = block.ghost;
+	block.nodex_end = block.nodex + block.ghost - 1;
+	Fluid1d* fluids = new Fluid1d[block.nx];
+
+	// then the interfaces reconstructioned value, N+1
+	Interface1d* interfaces = new Interface1d[block.nx + 1];
+	// then the flux, which the number is identical to interfaces
+	Flux1d** fluxes = Setflux_array(block);
+	//end the allocate memory part
+
+	//bulid or read mesh,
+	//since the mesh is all structured from left and right,
+	//there is no need for buliding the complex topology between cells and interfaces
+	//just using the index for address searching
+
+	block.left = 0.0;
+	block.right = 2.0;
+	block.dx = (block.right - block.left) / block.nodex;
+	//set the uniform geometry information
+	SetUniformMesh(block, fluids, interfaces, fluxes);
+
+
+	//初始化一下
+	ICfor_sinwave(fluids, block);
+
+	//初始化时间和时间步长
+	block.t = 0;//the current simulation time
+	block.step = 0; //the current step
+	int inputstep = 1;//输入一个运行的时间步数
+
+	runtime.finish_initial = clock(); //记录一下初始化需要的cpu时间....
+
+	while (block.t < tstop) //显式时间步循环，from t^n to t^n+1
+	{
+		//输入一个运行时间步数的主要好处是，你可以在开发阶段，先让程序跑几步停一下看看结果
+		//避免跑了一万步，发现从第一部就爆掉了
+		if (block.step % inputstep == 0)
+		{
+			cout << "pls cin interation step, if input is 0, then the program will exit " << endl;
+			cin >> inputstep;
+			if (inputstep == 0)
+			{
+				output1d(fluids, block); //这个是输出一维的流场数据
+				break;
+			}
+		}
+
+		//记录运算时间
+		if (runtime.start_compute == 0.0)
+		{
+			runtime.start_compute = clock();
+			cout << "runtime-start " << endl;
+		}
+		//end
+
+		//Copy the fluid vales to fluid old
+		CopyFluid_new_to_old(fluids, block);
+
+		//determine the cfl condtion
+		block.dt = Get_CFL(block, fluids, tstop);
+
+
+		for (int i = 0; i < block.stages; ++i)
+		{
+			//边界条件
+			leftboundary(fluids, block, bcvalue[0]);
+			rightboundary(fluids, block, bcvalue[1]);
+
+			//重构
+			Reconstruction_within_cell(interfaces, fluids, block);
+			Reconstruction_forg0(interfaces, fluids, block);
+
+			//计算通量
+			Calculate_flux(fluxes, interfaces, block, i);
+
+			//更新
+			Update(fluids, fluxes, block, i);
+		}
+
+		block.step++;
+		block.t = block.t + block.dt;
+		if ((block.t - tstop) > 0)
+		{
+			error_for_sinwave(fluids, block, tstop, error);
+		}
+	}
+	runtime.finish_compute = clock();
+	cout << "initializiation time is " << (float)(runtime.finish_initial - runtime.start_initial) / CLOCKS_PER_SEC << "seconds" << endl;
+	cout << "computational time is " << (float)(runtime.finish_compute - runtime.start_compute) / CLOCKS_PER_SEC << "seconds" << endl;
+	output1d(fluids, block);
+}
+
+void ICfor_sinwave(Fluid1d* fluids, Block1d block)
+{
+#pragma omp parallel for 
+	for (int i = block.nodex_begin; i <= block.nodex_end; i++)
+	{
+		double pi = 3.14159265358979323846;
+		fluids[i].primvar[0]
+			= 1.0 - 0.2 / pi / block.dx *
+			(cos(pi * (i + 1 - block.ghost) * block.dx) - cos(pi * (i - block.ghost) * block.dx));
+		fluids[i].primvar[1] = 1;
+		fluids[i].primvar[2] = 1;
+		Primvar_to_convar_1D(fluids[i].convar, fluids[i].primvar);
+
+	}
+}
+
+void error_for_sinwave(Fluid1d* fluids, Block1d block, double tstop, double* error)
+{
+
+	if (abs(block.t - tstop) <= (1e-10))
+	{
+
+		double error1 = 0;
+		double error2 = 0;
+		double error3 = 0;
+
+		for (int i = block.ghost; i < block.nx - block.ghost; i++)
+		{
+			double pi = 3.14159265358979323846;
+			int index = i;
+			double primvar0 = 1 - 0.2 / pi / block.dx * (cos(pi * (i + 1 - block.ghost) * block.dx) - cos(pi * (i - block.ghost) * block.dx));
+			error1 = error1 + abs(fluids[index].convar[0] - primvar0);
+			error2 = error2 + pow(fluids[index].convar[0] - primvar0, 2);
+			error3 = (error3 > abs(fluids[index].convar[0] - primvar0)) ? error3 : abs(fluids[index].convar[0] - primvar0);
+		}
+		error1 /= block.nodex;
+		error2 = sqrt(error2 / block.nodex);
+
+		error[0] = error1; error[1] = error2; error[2] = error3;
+
+		cout << scientific << "L1 error=" << error1 << " L2 error=" << error2 << " Linf error=" << error3 << endl;
+
+	}
+
+
+}
+
+void accuracy_sinwave_2d()
+{
+	int mesh_set = 3;
+	int mesh_number_start = 10;
+	double length = 2.0;
+	double CFL = 0.5;
+	double dt_ratio = 1.0;
+
+	double mesh_size_start = length / mesh_number_start;
+	int* mesh_number = new int[mesh_set];
+	double* mesh_size = new double[mesh_set];
+	double** error = new double* [mesh_set];
+
+	for (int i = 0; i < mesh_set; ++i)
+	{
+		error[i] = new double[3];
+		mesh_size[i] = mesh_size_start / pow(2, i);
+		mesh_number[i] = mesh_number_start * pow(2, i);
+		sinwave_2d(CFL, dt_ratio, mesh_number[i], error[i]);
+	}
+
+	output_error_form(CFL, dt_ratio, mesh_set, mesh_number, error);
+}
+
+void sinwave_2d(double& CFL, double& dt_ratio, int& mesh_number, double* error)
+{
+
+	Runtime runtime;
+	runtime.start_initial = clock();
+	Block2d block; // Class, geometry variables
+	block.uniform = false;
+	block.nodex = mesh_number;
+	block.nodey = mesh_number;
+	block.ghost = 3; // 5th-order reconstruction, should 3 ghost cell
+
+	double tstop = 2;
+	block.CFL = CFL;
+
+
+	K = 3; // 1d K=4, 2d K=3, 3d K=2;
+	Gamma = 1.4; // diatomic gas, r=1.4
+
+	//this part should rewritten ad gks2dsolver blabla
+	gks2dsolver = gks2nd_2d; // Emumeration, choose solver type
+	// 2nd means spatical 2nd-order, the traditional GKS; 2d means two dimensions
+
+	tau_type = Euler; // Emumeration, choose collision time tau type
+	// for accuracy test case, the tau type should be Euler and more accurately ZERO
+	// for inviscid case, use Euler type tau; for viscous case, use NS type tau
+
+	c1_euler = 0.0; // first coefficient in Euler type tau calculation
+	c2_euler = 0.0; // second coefficient in Euler type tau calculation
+	// when Euler type, tau is always zero;
+	// if the c1 c2 were given, (c1, c2, should be given together), the tau_num is determined by c1 c2, (c1*dt + c2*deltaP)
+	// if c1, c2 are both zero, (or not given value either), tau_num is also zero, which might be useful for accuracy test.
+	// when NS type, tau is always physical tau, relating to Mu or Nu;
+	// besides, tau_num would be determined by tau and c2 part, (tau + c2*deltaP),
+	// so, c2 should be given (c1 is useless now), and Mu or Nu should be given ONLY ONE;
+	// for NS type, specially, if Smooth is true, tau_num is determined ONLY by tau, (tau_num = tau), which means c2 is zero
+
+	Mu = 0.0;
+	Pr = 0.73;
+	R_gas = 1;
+
+	//prepare the boundary condtion function
+	Fluid2d* bcvalue = new Fluid2d[4]; // Class, primitive variables only for boundary
+	BoundaryCondition2dnew leftboundary(0); // 声明 边界条件的 函数指针
+	BoundaryCondition2dnew rightboundary(0); // 声明 边界条件的 函数指针
+	BoundaryCondition2dnew downboundary(0); // 声明 边界条件的 函数指针
+	BoundaryCondition2dnew upboundary(0); // 声明 边界条件的 函数指针
+
+	leftboundary = periodic_boundary_left;  // 给函数指针 赋值
+	rightboundary = periodic_boundary_right;  // 给函数指针 赋值
+	downboundary = periodic_boundary_down;  // 给函数指针 赋值
+	upboundary = periodic_boundary_up;  // 给函数指针 赋值
+
+	//prepare the reconstruction
+	gausspoint = 2; // fifth-order or sixth-order use THREE gauss points
+	// WENO5 has the function relating to arbitrary gausspoints
+	// WENO5_AO supports 2 gausspoint now, so fourth-order at most for spacial reconstruction (enough for two step fourth-order GKS)
+	SetGuassPoint(); // Function, set Gauss points coordinates and weight factor
+
+	reconstruction_variable = conservative; // Emumeration, choose the variables used for reconstruction type
+	wenotype = linear; // Emumeration, choose reconstruction type
+
+	cellreconstruction_2D_normal = WENO5_normal;  // reconstruction in normal directon
+	cellreconstruction_2D_tangent = WENO5_tangent;  // reconstruction in tangential directon
+	g0reconstruction_2D_normal = Center_5th_normal;  // reconstruction for g0 in normal directon
+	g0reconstruction_2D_tangent = Center_5th_tangent;  // reconstruction for g0 in tangential directon
+
+
+	//prepare the flux function
+	flux_function_2d = GKS2D_smooth; // 给函数指针 赋值 flux calculation type
+	// solver is GKS
+	// GKS2D means full solver (used for shock), GKS2D_smooth means smooth solver (used for non_shock)
+	// for accuracy test case, tau is ZERO, so GKS2D is equal to GKS2D_smooth for results, but the later is faster
+	// for inviscid flow, Euler type tau equals artifical viscosity;
+	// for viscous flow, NS type tau equals the modified collision time;
+	//end
+
+	//prepare time marching stratedgy
+	//time coe list must be 2d, end by _2D
+	timecoe_list_2d = S2O4_2D; // 两步四阶
+	Initial_stages(block);
+
+
+	// allocate memory for 2-D fluid field
+	// in a standard finite element method, we have
+	// first the cell average value, N*N
+	block.nx = block.nodex + 2 * block.ghost;
+	block.ny = block.nodey + 2 * block.ghost;
+
+	Fluid2d* fluids = Setfluid(block); // Function, input a class (geometry), output the pointer of one class (conservative variables)
+	// then the interfaces reconstructioned value, (N+1)*(N+1)
+	Interface2d* xinterfaces = Setinterface_array(block);
+	Interface2d* yinterfaces = Setinterface_array(block);
+	// then the flux, which the number is identical to interfaces
+	Flux2d_gauss** xfluxes = Setflux_gauss_array(block);
+	Flux2d_gauss** yfluxes = Setflux_gauss_array(block);
+	//end the allocate memory part
+
+	block.left = 0.0;
+	block.right = 2.0;
+	block.down = 0.0;
+	block.up = 2.0;
+	block.dx = (block.right - block.left) / block.nodex;
+	block.dy = (block.up - block.down) / block.nodey;
+	block.overdx = 1 / block.dx;
+	block.overdy = 1 / block.dy;
+	//set the uniform geometry information
+	SetUniformMesh(block, fluids, xinterfaces, yinterfaces, xfluxes, yfluxes); // Function, set mesh
+
+	//end
+
+	ICfor_sinwave_2d(fluids, block);
+
+	runtime.finish_initial = clock();
+	block.t = 0; //the current simulation time
+	block.step = 0; //the current step
+
+	int inputstep = 1;//input a certain step
+
+	while (block.t < tstop)
+	{
+
+		if (block.step % inputstep == 0)
+		{
+			cout << "pls cin interation step, if input is 0, then the program will exit " << endl;
+			cin >> inputstep;
+			if (inputstep == 0)
+			{
+				output2d_binary(fluids, block); // 二进制输出，用单元中心点代替单元
+				break;
+			}
+		}
+		if (runtime.start_compute == 0.0)
+		{
+			runtime.start_compute = clock();
+			cout << "runtime-start " << endl;
+		}
+
+		CopyFluid_new_to_old(fluids, block); // Function, copy variables
+		//determine the cfl condtion
+		block.dt = Get_CFL(block, fluids, tstop); // Function, get real CFL number
+
+		if ((block.t + block.dt - tstop) > 0)
+		{
+			block.dt = tstop - block.t + 1e-16;
+		}
+
+		for (int i = 0; i < block.stages; i++)
+		{
+			//after determine the cfl condition, let's implement boundary condtion
+			leftboundary(fluids, block, bcvalue[0]);
+			rightboundary(fluids, block, bcvalue[1]);
+			downboundary(fluids, block, bcvalue[2]);
+			upboundary(fluids, block, bcvalue[3]);
+			//cout << "after bc " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+
+			Convar_to_Primvar(fluids, block); // Function
+			//then is reconstruction part, which we separate the left or right reconstrction
+			//and the center reconstruction
+			Reconstruction_within_cell(xinterfaces, yinterfaces, fluids, block); // Function
+			//cout << "after cell recon " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+
+			Reconstruction_forg0(xinterfaces, yinterfaces, fluids, block); // Function
+			//cout << "after g0 recon " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+
+			//then is solver part
+			Calculate_flux(xfluxes, yfluxes, xinterfaces, yinterfaces, block, i); // Function
+			//cout << "after flux calcu " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+
+			//then is update flux part
+			Update(fluids, xfluxes, yfluxes, block, i); // Function
+			//cout << "after updating " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+
+		}
+
+		block.step++;
+		block.t = block.t + block.dt;
+
+		if (block.step % 1000 == 0)
+		{
+			cout << "step 1000 time is " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+		}
+		if ((block.t - tstop) >= 0)
+		{
+			output2d_binary(fluids, block);
+		}
+
+	}
+
+	error_for_sinwave_2d(fluids, block, tstop, error); // Function, get error
+	runtime.finish_compute = clock();
+
+	cout << "the total run time is " << (double)(runtime.finish_compute - runtime.start_initial) / CLOCKS_PER_SEC << " second !" << endl;
+	cout << "initializing time is" << (double)(runtime.finish_initial - runtime.start_initial) / CLOCKS_PER_SEC << " second !" << endl;
+	cout << "the pure computational time is" << (double)(runtime.finish_compute - runtime.start_compute) / CLOCKS_PER_SEC << " second !" << endl;
+
+}
+
+void ICfor_sinwave_2d(Fluid2d* fluids, Block2d block)
+{
+	for (int i = block.ghost; i < block.nx - block.ghost; i++)
+	{
+		for (int j = block.ghost; j < block.ny - block.ghost; j++)
+		{
+			double pi = 3.14159265358979323846;
+			int index = i * (block.ny) + j;
+			double xleft = (i - block.ghost) * block.dx;
+			double xright = (i + 1 - block.ghost) * block.dx;
+			double yleft = (j - block.ghost) * block.dy;
+			double yright = (j + 1 - block.ghost) * block.dy;
+
+			//case in two dimensional x-y-plane
+			double k1 = sin(pi * (xright + yright));
+			double k2 = sin(pi * (xright + yleft));
+			double k3 = sin(pi * (xleft + yright));
+			double k4 = sin(pi * (xleft + yleft));
+			fluids[index].primvar[0] = 1.0 - 0.2 / pi / pi / block.dx / block.dy * ((k1 - k2) - (k3 - k4));
+
+			fluids[index].primvar[1] = 1;
+			fluids[index].primvar[2] = 1;
+			fluids[index].primvar[3] = 1;
+		}
+	}
+#pragma omp parallel for
+	for (int i = block.ghost; i < block.nx - block.ghost; i++)
+	{
+		for (int j = block.ghost; j < block.ny - block.ghost; j++)
+		{
+			int index = i * (block.ny) + j;
+			Primvar_to_convar_2D(fluids[index].convar, fluids[index].primvar);
+		}
+	}
+}
+
+void error_for_sinwave_2d(Fluid2d* fluids, Block2d block, double tstop, double* error)
+{
+	if (abs(block.t - tstop) <= (1e-10))
+	{
+		cout << "Accuracy-residual-file-output" << endl;
+		double error1 = 0;
+		double error2 = 0;
+		double error3 = 0;
+		for (int i = block.ghost; i < block.nx - block.ghost; i++)
+		{
+			for (int j = block.ghost; j < block.ny - block.ghost; j++)
+			{
+				double pi = 3.14159265358979323846;
+				int index = i * block.ny + j;
+				double xleft = (i - block.ghost) * block.dx;
+				double xright = (i + 1 - block.ghost) * block.dx;
+				double yleft = (j - block.ghost) * block.dy;
+				double yright = (j + 1 - block.ghost) * block.dy;
+
+				//case in two dimensional x-y-plane
+				double k1 = sin(pi * (xright + yright));
+				double k2 = sin(pi * (xright + yleft));
+				double k3 = sin(pi * (xleft + yright));
+				double k4 = sin(pi * (xleft + yleft));
+				double primvar0 = 1.0 - 0.2 / pi / pi / block.dx / block.dy * ((k1 - k2) - (k3 - k4));
+
+
+				error1 = error1 + abs(fluids[index].convar[0] - primvar0);
+				error2 = error2 + pow(fluids[index].convar[0] - primvar0, 2);
+				error3 = (error3 > abs(fluids[index].convar[0] - primvar0)) ? error3 : abs(fluids[index].convar[0] - primvar0);
+			}
+		}
+		error1 /= (block.nodex * block.nodey);
+		error2 = sqrt(error2 / (block.nodex * block.nodey));
+		error[0] = error1; error[1] = error2; error[2] = error3;
+	}
+}

code/accuracy_test.h

@@ -0,0 +1,11 @@
+#pragma once
+#include"fvm_gks2d.h"
+#include"output.h"
+void accuracy_sinwave_1d();
+void accuracy_sinwave_1d(double& CFL, double& dt_ratio, int& mesh_number, double* error);
+void ICfor_sinwave(Fluid1d* fluids, Block1d block);
+void error_for_sinwave(Fluid1d* fluids, Block1d block, double tstop, double* error);
+void accuracy_sinwave_2d();
+void sinwave_2d(double &CFL, double &dt_ratio, int& mesh_number, double *error);
+void ICfor_sinwave_2d(Fluid2d *fluids, Block2d block);
+void error_for_sinwave_2d(Fluid2d *fluids, Block2d block, double tstop,double *error);

code/boundary_layer.cpp

@@ -0,0 +1,373 @@
+#include"boundary_layer.h"
+
+
+void boundary_layer()
+{
+
+	Runtime runtime;
+	runtime.start_initial = clock();
+	Block2d block;
+	block.uniform = false;
+	block.ghost = 3;
+
+	double tstop = 20000;
+	block.CFL = 0.3;
+
+	K = 3;
+	Gamma = 1.4;
+	Pr = 1.0;
+
+	double renum = 1e5;
+	double den_ref = 1.0;
+	double u_ref = 0.15;
+	double L = 100;
+
+	//prepare the boundary condtion function
+	Fluid2d* bcvalue = new Fluid2d[1];
+	bcvalue[0].primvar[0] = den_ref;
+	bcvalue[0].primvar[1] = u_ref;
+	bcvalue[0].primvar[2] = 0.0;
+	bcvalue[0].primvar[3] = den_ref / Gamma;
+
+	Mu = den_ref * u_ref * L / renum;
+	cout << Mu << endl;
+
+
+
+	//prepare the reconstruction function
+
+	gausspoint = 1;
+	SetGuassPoint();
+	reconstruction_variable = conservative;
+	wenotype = wenoz;//只有WENO才会起作用
+	cellreconstruction_2D_normal = Vanleer_normal;
+	cellreconstruction_2D_tangent = Vanleer_tangent;
+	g0reconstruction_2D_normal = Center_3rd_normal;
+	g0reconstruction_2D_tangent = Center_3rd_tangent;
+
+	is_reduce_order_warning = true;
+
+	flux_function_2d = GKS2D;
+	gks2dsolver = gks2nd_2d;
+	tau_type = NS;
+	c1_euler = 0.0;
+	c2_euler = 0;
+
+	timecoe_list_2d = S1O2_2D;
+
+
+	Initial_stages(block);
+
+
+	Fluid2d* fluids = NULL;
+	Interface2d* xinterfaces = NULL;
+	Interface2d* yinterfaces = NULL;
+	Flux2d_gauss** xfluxes = NULL;
+	Flux2d_gauss** yfluxes = NULL;
+
+	//读入2维结构化网格，先找到网格名字
+	string grid = add_mesh_directory_modify_for_linux()
+		+ "structured-mesh/flat-plate-str.plt";
+
+	//读入2维结构化网格，格式是tecplot的3维binary，z方向只有一层。可以通过pointwise输出tecplot的binary格式实现
+	Read_structured_mesh
+	(grid, &fluids, &xinterfaces, &yinterfaces, &xfluxes, &yfluxes, block);
+
+	ICfor_uniform_2d(fluids, bcvalue[0].primvar, block);
+
+	runtime.finish_initial = clock();
+	block.t = 0;
+	block.step = 0;
+	int inputstep = 1;
+
+	while (block.t < tstop)
+	{
+
+		if (block.step % inputstep == 0)
+		{
+			cout << "pls cin interation step, if input is 0, then the program will exit " << endl;
+			cin >> inputstep;
+			if (inputstep == 0)
+			{
+				output2d(fluids, block);
+				break;
+			}
+		}
+		if (runtime.start_compute == 0.0)
+		{
+			runtime.start_compute = clock();
+			cout << "runtime-start " << endl;
+		}
+
+		CopyFluid_new_to_old(fluids, block);
+		block.dt = Get_CFL(block, fluids, tstop);
+
+		for (int i = 0; i < block.stages; i++)
+		{
+
+			boundaryforBoundary_layer(fluids, block, bcvalue[0]);
+
+			Convar_to_Primvar(fluids, block);
+
+			Reconstruction_within_cell(xinterfaces, yinterfaces, fluids, block);
+
+			Reconstruction_forg0(xinterfaces, yinterfaces, fluids, block);
+
+			Calculate_flux(xfluxes, yfluxes, xinterfaces, yinterfaces, block, i);
+
+			Update(fluids, xfluxes, yfluxes, block, i);
+
+		}
+
+		++block.step;
+		block.t = block.t + block.dt;
+
+		if (block.step % 100 == 0)
+		{
+			cout << "step 10 time is " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+		}
+		Residual2d(fluids, block, 10);
+		output_wall_info_boundary_layer(xinterfaces, yinterfaces, xfluxes, yfluxes, block, L);
+		output_skin_friction_boundary_layer(fluids, yinterfaces, yfluxes, block, u_ref, renum, L);
+		if (block.step % 5000 == 0 || abs(block.t - tstop) < 1e-14)
+		{
+			output2d(fluids, block);
+			char loc1[] = "0050";
+			char loc2[] = "0105";
+			char loc3[] = "0165";
+			char loc4[] = "0265";
+			char loc5[] = "0305";
+			char loc6[] = "0365";
+			char loc7[] = "0465";
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.050 * L, loc1);
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.105 * L, loc2);
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.165 * L, loc3);
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.265 * L, loc4);
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.305 * L, loc5);
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.365 * L, loc6);
+			output_boundary_layer(fluids, block, bcvalue[0].primvar, 0.465 * L, loc7);
+
+		}
+
+	}
+	runtime.finish_compute = clock();
+	;
+	cout << "\n the total running time is " << (double)(runtime.finish_compute - runtime.start_initial) / CLOCKS_PER_SEC << "秒！" << endl;
+	cout << "\n the time for initializing is " << (double)(runtime.finish_initial - runtime.start_initial) / CLOCKS_PER_SEC << "秒！" << endl;
+	cout << "\n the time for computing is " << (double)(runtime.finish_compute - runtime.start_compute) / CLOCKS_PER_SEC << "秒！" << endl;
+
+}
+
+void boundaryforBoundary_layer(Fluid2d* fluids, Block2d block, Fluid2d bcvalue)
+{
+
+	inflow_boundary_left(fluids, block, bcvalue);
+
+	inflow_boundary_right(fluids, block, bcvalue);
+	//free_boundary_up(fluids, block, bcvalue);
+	//outflow_boundary_up(fluids, block, bcvalue);
+	inflow_boundary_up(fluids, block, bcvalue);
+
+	//the boundary at the bottom
+	int order = block.ghost;
+
+	for (int j = order - 1; j >= 0; j--)
+	{
+		for (int i = order; i < block.nx; i++)
+		{
+			int index = i * block.ny + j;
+			int ref = i * block.ny + 2 * order - 1 - j;
+			if (fluids[index].coordx < 0.0)
+			{
+
+				fluids[index].convar[0] = fluids[ref].convar[0];
+				fluids[index].convar[1] = fluids[ref].convar[1];
+				fluids[index].convar[2] = -fluids[ref].convar[2];
+				fluids[index].convar[3] = fluids[ref].convar[3];
+
+
+
+			}
+			else
+			{
+				fluids[index].convar[0] = fluids[ref].convar[0];
+				fluids[index].convar[1] = -fluids[ref].convar[1];
+				fluids[index].convar[2] = -fluids[ref].convar[2];
+				fluids[index].convar[3] = fluids[ref].convar[3];
+
+
+			}
+
+
+		}
+
+	}
+}
+
+void output_boundary_layer(Fluid2d* fluids, Block2d block, double* bcvalue, double coordx, char* label)
+{
+	int order = block.ghost;
+	int i;
+	for (int k = 0; k < block.nx; k++)
+	{
+		int index = k * block.ny + order;
+		if (fluids[index].coordx > coordx)
+		{
+			i = k;
+			break;
+		}
+	}
+
+	stringstream name;
+	name << "result/boundarylayer_line" << label << setfill('0') << setw(5) << "step" << block.step << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << "variables =ys,density,us,vs,pressure" << endl;
+	out << "zone i = " << 1 << ",j = " << block.nodey << ", F=POINT" << endl;
+
+	double vel = bcvalue[1];
+	double nu = Mu / bcvalue[0];
+	for (int j = order; j < block.ny - order; j++)
+	{
+		int index = i * block.ny + j;
+		Convar_to_primvar_2D(fluids[index].primvar, fluids[index].convar);
+
+		out << fluids[index].coordy / (sqrt(nu * fluids[index].coordx / vel)) << " "
+			<< fluids[index].primvar[0] << " "
+			<< fluids[index].primvar[1] / vel << " "
+			<< fluids[index].primvar[2] / (sqrt(nu * vel / fluids[index].coordx)) << " "
+			<< fluids[index].primvar[3] << " "
+			<< endl;
+	}
+
+	out.close();
+
+}
+
+void output_wall_info_boundary_layer(Interface2d* xInterfaces, Interface2d* yInterfaces,
+	Flux2d_gauss** xFluxes, Flux2d_gauss** yFluxes, Block2d block, double L_ref)
+{
+	//there are four blocks
+	int outstep = 20;
+	if (block.step % outstep != 0)
+	{
+		return;
+	}
+	double t = block.t;
+	double dt = block.timecoefficient[block.stages - 1][0][0] * block.dt;
+	double var[8];
+	Array_zero(var, 8);
+
+
+	//the contribute from the bottom
+	for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+	{
+
+		int ref = i * (block.ny + 1) + block.ghost;
+
+		if (yInterfaces[ref].x > 0.0)
+		{
+			double f_local[4];
+			double derf_local[4];
+			for (int m = 0; m < gausspoint; m++)
+			{
+				Local_to_Global(f_local, yFluxes[ref][0].gauss[m].f, yFluxes[ref][0].gauss[m].normal);
+				Local_to_Global(derf_local, yFluxes[ref][0].gauss[m].derf, yFluxes[ref][0].gauss[m].normal);
+			}
+			double l = yInterfaces[ref].length;
+			double sign = 1;
+			for (int k = 0; k < 4; k++)
+			{
+				for (int m = 0; m < gausspoint; m++)
+				{
+					var[k] += sign * l * yFluxes[ref][block.stages - 1].gauss[m].x[k] / dt;
+					var[k + 4] += sign * l * gauss_weight[m] * (f_local[k]) / block.dt;
+				}
+			}
+		}
+
+	}
+
+	ofstream out;
+	if (block.step == outstep)
+	{
+		out.open("result/boundary-layer-force-history.plt", ios::out);
+	}
+	else
+	{
+		out.open("result/boundary-layer-force-history.plt", ios::ate | ios::out | ios::in);
+	}
+
+	if (!out.is_open())
+	{
+		cout << "cannot find square-cylinder-force-history.plt" << endl;
+		cout << "a new case will start" << endl;
+	}
+
+	if (block.step == outstep)
+	{
+		out << "# CFL number is " << block.CFL << endl;
+		out << "variables = t,mass_transfer,cd,cl,heatflux,";
+		out << "mass_transfer_instant, cd_instant, cl_instant, heatflux_instant" << endl;
+	}
+
+	Array_scale(var, 8, L_ref);
+
+	out << block.t << " "
+		<< var[0] << " " << var[1] << " "
+		<< var[2] << " " << var[3] << " ";
+	out << var[4] << " " << var[5] << " "
+		<< var[6] << " " << var[7] << " ";
+	out << endl;
+	out.close();
+}
+
+
+void output_skin_friction_boundary_layer
+(Fluid2d* fluids, Interface2d* yInterfaces, Flux2d_gauss** yFluxes, Block2d block, double u_ref, double re_ref, double L_ref)
+{
+	int outstep = 1000;
+	if (block.step % outstep != 0)
+	{
+		return;
+	}
+	double t = block.t;
+	double dt = block.timecoefficient[0][0][0] * block.dt;
+	double var[3];
+	Array_zero(var, 3);
+
+	ofstream out;
+	out.open("result/boundary-layer-skin-friction.plt", ios::out);
+	out << "variables = log<sub>10</sub>Re<sub>x</sub>, log<sub>10</sub>C<sub>f</sub>, C<sub>p</sub>";
+	out << endl;
+	//output the data
+	//the contribute from the bottom
+	for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+	{
+		int ref = i * (block.ny + 1) + block.ghost;
+		int ref_cell = i * (block.ny) + block.ghost;
+		if (yInterfaces[ref].x > 0.0)
+		{
+			Array_zero(var, 3);
+			for (int k = 1; k < 3; k++)
+			{
+				for (int m = 0; m < gausspoint; m++)
+				{
+					var[k] += yFluxes[ref][0].gauss[m].x[k] / dt;
+				}
+			}
+			var[1] = -var[1];
+			var[2] = (var[2] - 1.0 / Gamma);
+			double rex = yInterfaces[ref].x / L_ref * re_ref;
+
+			for (int k = 0; k < 3; ++k)
+			{
+				var[k] /= (0.5 * u_ref * u_ref);
+			}
+			out << log10(rex) << " " << log10(var[1]) << " " << var[2] << endl;
+
+		}
+	}
+	out.close();
+}

code/boundary_layer.h

@@ -0,0 +1,11 @@
+#pragma once
+#include"fvm_gks2d.h"
+#include"output.h"
+#include"mesh_read.h"
+void boundary_layer();
+void boundaryforBoundary_layer(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void output_boundary_layer(Fluid2d* fluids, Block2d block, double* bcvalue, double coordx, char* label);
+void output_wall_info_boundary_layer(Interface2d* xInterfaces, Interface2d* yInterfaces,
+	Flux2d_gauss** xFluxes, Flux2d_gauss** yFluxes, Block2d block, double L_ref);
+void output_skin_friction_boundary_layer
+(Fluid2d* fluids, Interface2d* yInterfaces, Flux2d_gauss** yFluxes, Block2d block, double u_ref, double re_ref, double L_ref);

code/cylinder.cpp

@@ -0,0 +1,148 @@
+#include"cylinder.h"
+
+void cylinder()
+{
+	Runtime runtime;
+	runtime.start_initial = clock();
+	Block2d block;
+	block.uniform = false;
+	block.ghost = 3;
+
+	double tstop = 1;
+	block.CFL = 0.5;
+
+	K = 3;
+	Gamma = 1.4;
+	Pr = 1.0;
+
+	double renum = 1e3;
+	double den_ref = 1.0;
+	double u_ref = 5.0;
+	double L = 1;  //cylinder diameter 1
+
+	Fluid2d icvalue;
+	icvalue.primvar[0] = den_ref;
+	icvalue.primvar[1] = u_ref;
+	icvalue.primvar[2] = 0.0;
+	icvalue.primvar[3] = den_ref / Gamma;
+
+	Mu = den_ref * u_ref * L / renum;
+	cout << Mu << endl;
+
+	gausspoint = 1;
+	SetGuassPoint();
+	reconstruction_variable = conservative;
+	wenotype = wenoz;
+
+
+	cellreconstruction_2D_normal = Vanleer_normal;
+	cellreconstruction_2D_tangent = Vanleer_tangent;
+	g0reconstruction_2D_normal = Center_do_nothing_normal;
+	g0reconstruction_2D_tangent = Center_all_collision_multi;
+
+	is_reduce_order_warning = false;
+	flux_function_2d = GKS2D;
+	gks2dsolver = gks2nd_2d;
+	tau_type = NS;
+	c1_euler = 0.05;
+	c2_euler = 5; //强间断 c2建议取大于等于5
+
+	timecoe_list_2d = S1O2_2D;
+
+	Initial_stages(block);
+
+	Fluid2d* fluids = NULL;
+	Interface2d* xinterfaces = NULL;
+	Interface2d* yinterfaces = NULL;
+	Flux2d_gauss** xfluxes = NULL;
+	Flux2d_gauss** yfluxes = NULL;
+
+	string grid = add_mesh_directory_modify_for_linux()
+		+ "structured-mesh/half-cylinder-str-1.plt";
+
+	Read_structured_mesh
+	(grid, &fluids, &xinterfaces, &yinterfaces, &xfluxes, &yfluxes, block);
+
+	ICfor_uniform_2d(fluids, icvalue.primvar, block);
+
+	runtime.finish_initial = clock();
+	block.t = 0;
+	block.step = 0;
+
+	int inputstep = 1;
+	while (block.t < tstop)
+	{
+		if (block.step % inputstep == 0)
+		{
+			cout << "pls cin interation step, if input is 0, then the program will exit " << endl;
+			cin >> inputstep;
+			if (inputstep == 0)
+			{
+				output2d_center(fluids, block);
+				break;
+			}
+		}
+		if (runtime.start_compute == 0.0)
+		{
+			runtime.start_compute = clock();
+			cout << "runtime-start " << endl;
+		}
+
+		CopyFluid_new_to_old(fluids, block);
+		Convar_to_Primvar(fluids, block);
+
+		block.dt = Get_CFL(block, fluids, tstop);
+
+		for (int i = 0; i < block.stages; i++)
+		{
+
+			boundary_for_cylinder(fluids, xinterfaces, block, icvalue);
+
+			Reconstruction_within_cell(xinterfaces, yinterfaces, fluids, block);
+
+			Reconstruction_forg0(xinterfaces, yinterfaces, fluids, block);
+
+			Calculate_flux(xfluxes, yfluxes, xinterfaces, yinterfaces, block, i);
+
+			Update(fluids, xfluxes, yfluxes, block, i);
+
+		}
+
+
+		++block.step;
+		block.t = block.t + block.dt;
+
+		if (block.step % 100 == 0)
+		{
+			cout << "step 10 time is " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+		}
+		Residual2d(fluids, block, 10);
+
+	}
+	runtime.finish_compute = clock();
+	cout << "\n the total running time is " << (double)(runtime.finish_compute - runtime.start_initial) / CLOCKS_PER_SEC << " second!" << endl;
+	cout << "\n the time for initializing is " << (double)(runtime.finish_initial - runtime.start_initial) / CLOCKS_PER_SEC << " second!" << endl;
+	cout << "\n the time for computing is " << (double)(runtime.finish_compute - runtime.start_compute) / CLOCKS_PER_SEC << " second!" << endl;
+
+	//output_prim_variable_at_final_time
+	output2d_center(fluids, block);
+}
+
+
+void boundary_for_cylinder(Fluid2d* fluids, Interface2d* faces, Block2d block, Fluid2d bcvalue)
+{
+	Fluid2d wall_vel;
+	wall_vel.primvar[0] = bcvalue.primvar[0];
+	wall_vel.primvar[1] = 0.0;
+	wall_vel.primvar[2] = 0.0;
+	wall_vel.primvar[3] = bcvalue.primvar[3];
+
+	noslip_adiabatic_boundary_right(fluids, faces, block, wall_vel);
+
+	inflow_boundary_left(fluids, block, bcvalue);
+
+	free_boundary_down(fluids, block, bcvalue);
+	free_boundary_up(fluids, block, bcvalue);
+
+
+}

code/cylinder.h

@@ -0,0 +1,6 @@
+#pragma once
+#include"fvm_gks2d.h"
+#include"output.h"
+#include"mesh_read.h"
+void cylinder();
+void boundary_for_cylinder(Fluid2d* fluids, Interface2d* faces, Block2d block, Fluid2d bcvalue);

code/function.cpp

@@ -0,0 +1,303 @@
+#include"function.h"
+
+int K = 0; //initialization  //双原子，r=1.4，K=（4-2r）/（r-1）=3
+double Gamma = 2.0 / (K + 2) + 1.0; //initialization
+
+double T_inf = 285; //the real temperature, 20 Celcius degree
+double R_gas = 8.31441 / (28.959e-3); //initialization //the gas constant for real ideal air
+double Mu = -1.0; // mu = rho*u*L/Re
+double Nu = -1.0; // Nu = u*L/Re
+TAU_TYPE tau_type=Euler; //无粘还是粘性问题
+//tau_type, euler==0,ns==1,sutherland==2, Power_law==3
+
+void Copy_Array(double* target, double* origin, int dim)
+{
+	for (int i = 0; i < dim; i++)
+	{
+		target[i] = origin[i];
+	}
+}
+void Array_zero(double* target, int dim)
+{
+	for (int i = 0; i < dim; i++)
+	{
+		target[i] = 0.0;
+	}
+}
+void Array_subtract(double* target, double* subtractor, int dim)
+{
+	for (int i = 0; i < dim; i++)
+	{
+		target[i] -= subtractor[i];
+	}
+}
+
+void Array_scale(double* target, int dim, double value)
+{
+	for (int i = 0; i < dim; i++)
+	{
+		target[i] *= value;
+	}
+}
+
+//从守恒变量计算原始变量
+void Convar_to_primvar_1D(double* primvar, double convar[3])
+{
+	primvar[0] = convar[0];
+	primvar[1] = U(convar[0], convar[1]);
+	primvar[2] = Pressure(convar[0], convar[1], convar[2]);
+
+}
+double Pressure(double density, double densityu, double densityE)
+{
+	return (Gamma - 1) * (densityE - 0.5 * densityu * densityu / density);
+}
+double Temperature(double density, double pressure)
+{
+	return pressure / R_gas / density;
+}
+double entropy(double density, double pressure)
+{
+	return log(pressure / pow(density, Gamma));
+}
+double Soundspeed(double density, double pressure)
+{
+	return sqrt(Gamma * pressure / density);
+}
+void Primvar_to_convar_1D(double* convar, double primvar[3])
+{
+	convar[0] = primvar[0];
+	convar[1] = DensityU(primvar[0], primvar[1]);
+	convar[2] = DensityE(primvar[0], primvar[1], primvar[2]);
+}
+double DensityU(double density, double u)
+{
+	return density * u;
+}
+double DensityE(double density, double u, double pressure)
+{
+	return density * (pressure / density / (Gamma - 1) + 0.5 * (u * u));
+}
+void Convar_to_char1D(double* character, double primvar[3], double convar[3])
+{
+	double c = sqrt(Gamma * primvar[2] / primvar[0]);
+	double	alfa = (Gamma - 1.0) / (2.0 * c * c);
+	double u = primvar[1];
+	double s[3][3];
+	s[0][0] = alfa * (0.5 * u * u + u * c / (Gamma - 1.0));
+	s[0][1] = alfa * (-u - c / (Gamma - 1.0));
+	s[0][2] = alfa;
+	s[1][0] = alfa * (-u * u + 2.0 * c * c / (Gamma - 1.0));
+	s[1][1] = alfa * 2.0 * u;
+	s[1][2] = -2.0 * alfa;
+	s[2][0] = alfa * (0.5 * u * u - u * c / (Gamma - 1.0));
+	s[2][1] = alfa * (-u + c / (Gamma - 1.0));
+	s[2][2] = alfa;
+
+	for (int i = 0; i < 3; i++)
+	{
+		character[i] = 0;
+	}
+	for (int i = 0; i < 3; i++)
+	{
+		for (int j = 0; j < 3; j++)
+		{
+			character[i] = character[i] + s[i][j] * convar[j];
+		}
+	}
+}
+void Char_to_convar1D(double* convar, double primvar[3], double charvar[3])
+{
+	double c = sqrt(Gamma * primvar[2] / primvar[0]);
+	double u = primvar[1];
+	double	h = 0.5 * u * u + c * c / (Gamma - 1.0);
+	double s[3][3];
+	s[0][0] = 1.0;
+	s[0][1] = 1.0;
+	s[0][2] = 1.0;
+	s[1][0] = u - c;
+	s[1][1] = u;
+	s[1][2] = u + c;
+	s[2][0] = h - u * c;
+	s[2][1] = u * u / 2.0;
+	s[2][2] = h + u * c;
+
+	for (int i = 0; i < 3; i++)
+	{
+		convar[i] = 0;
+		for (int j = 0; j < 3; j++)
+		{
+			convar[i] = convar[i] + s[i][j] * charvar[j];
+		}
+	}
+}
+
+double U(double density, double q_densityu)
+{
+	return q_densityu / density;
+}
+double V(double density, double q_densityv)
+{
+	return q_densityv / density;
+}
+double Pressure(double density, double q_densityu, double q_densityv, double q_densityE)
+{
+	return (Gamma - 1.0)*(q_densityE - 0.5*q_densityu*q_densityu / density - 0.5*q_densityv*q_densityv / density);
+}
+
+
+
+
+//从原始变量计算守恒变量
+double Q_densityu(double density, double u)
+{
+	double q_densityu = density*u;
+	return q_densityu;
+}
+double Q_densityv(double density, double v)
+{
+	double q_densityv = density*v;
+	return q_densityv;
+}
+double Q_densityE(double density, double u, double v, double pressure)
+{
+	double q_densityE = density*(pressure / density / (Gamma - 1) + 0.5*(u*u + v*v));
+	return q_densityE;
+}
+
+
+
+void Primvar_to_convar_2D(double *convar, double primvar[4])//2d version
+{
+    convar[0] = primvar[0];
+    convar[1] = Q_densityu(primvar[0],primvar[1]);
+    convar[2] = Q_densityv(primvar[0],primvar[2]);
+    convar[3] = Q_densityE(primvar[0], primvar[1], primvar[2], primvar[3]);
+}
+void Convar_to_primvar_2D(double *primvar, double convar[4])//2d version //density, U, V, pressure
+{
+    primvar[0] = convar[0];
+	primvar[1] = U(convar[0], convar[1]);
+	primvar[2] = V(convar[0], convar[2]);
+	primvar[3] = Pressure(convar[0], convar[1], convar[2], convar[3]);
+}
+
+void Convar_to_char(double *character, double *primvar, double convar[4])
+{
+	double c = sqrt(Gamma *primvar[3] / primvar[0]);
+	double overc = 1.0 / c;
+	double	alfa = (Gamma - 1.0)*0.5*overc*overc;
+	double u = primvar[1];
+	double v = primvar[2];
+	double s[4][4];
+	s[0][0] = 1.0 - alfa*(u*u + v*v);
+	s[0][1] = 2.0*alfa*u;
+	s[0][2] = 2.0*alfa*v;
+	s[0][3] = -2.0*alfa;
+	s[1][0] = -v;
+	s[1][1] = 0.;
+	s[1][2] = 1.;
+	s[1][3] = 0.;
+	s[2][0] = 0.5*(-u *overc + alfa*(u*u + v*v));
+	s[2][1] = 0.5 *overc - alfa*u;
+	s[2][2] = -alfa*v;
+	s[2][3] = alfa;
+	s[3][0] = 0.5*(u *overc + alfa*(u*u + v*v));
+	s[3][1] = -0.5 *overc - alfa*u;
+	s[3][2] = -alfa*v;
+	s[3][3] = alfa;
+	for (int i = 0; i < 4; i++)
+	{
+		character[i] = 0;
+	}
+	for (int i = 0; i < 4; i++)
+	{
+		for (int j = 0; j < 4; j++)
+		{
+			character[i] = character[i] + s[i][j] * convar[j];
+		}
+	}
+}
+void Char_to_convar(double *convar, double *primvar, double character[4])
+{
+	double c = sqrt(Gamma *primvar[3] / primvar[0]);
+	double	alfa = (Gamma - 1.0) / (2.0*c*c);
+	double u = primvar[1];
+	double v = primvar[2];
+	double s[4][4];
+	s[0][0] = 1.;
+	s[0][1] = 0.0;
+	s[0][2] = 1.0;
+	s[0][3] = 1.0;
+	s[1][0] = u;
+	s[1][1] = 0.;
+	s[1][2] = u + c;
+	s[1][3] = u - c;
+	s[2][0] = v;
+	s[2][1] = 1.;
+	s[2][2] = v;
+	s[2][3] = v;
+	s[3][0] = 0.5*(u*u + v*v);
+	s[3][1] = v;
+	s[3][2] = 0.5*(u*u + v*v) + c*u + c*c / (Gamma - 1.0);
+	s[3][3] = 0.5*(u*u + v*v) - c*u + c*c / (Gamma - 1.0);
+
+	for (int i = 0; i < 4; i++)
+	{
+		convar[i] = 0;
+		for (int j = 0; j < 4; j++)
+		{
+			convar[i] = convar[i] + s[i][j] * character[j];
+		}
+	}
+}
+
+void Global_to_Local(double* change, double *normal)
+{
+	double temp[2];
+	temp[0] = change[1];
+	temp[1] = change[2];
+	change[1] = temp[0] * normal[0] + temp[1] * normal[1];
+	change[2] = -temp[0] * normal[1] + temp[1] * normal[0];
+}
+
+void Global_to_Local(double* change, double* origin, double* normal)
+{
+	change[0] = origin[0];
+	change[1] = origin[1] * normal[0] + origin[2] * normal[1];
+	change[2] = -origin[1] * normal[1] + origin[2] * normal[0];
+	change[3] = origin[3];
+}
+
+void Local_to_Global(double *change,double *normal)
+{
+    double temp[2];
+    temp[0] = change[1];
+    temp[1] = change[2];
+    change[1] = temp[0] * normal[0] - temp[1] * normal[1];
+    change[2] = temp[0] * normal[1] + temp[1] * normal[0];
+}
+
+void Local_to_Global(double* change, double* origin, double *normal)
+{
+	change[0] = origin[0];
+	change[1] = origin[1] * normal[0] - origin[2] * normal[1];
+	change[2] = origin[1] * normal[1] + origin[2] * normal[0];
+	change[3] = origin[3];
+}
+
+double Dtx(double dtx, double dx, double CFL, double density, double u,double v, double pressure)
+{
+	double tmp;
+	tmp = sqrt(u*u+v*v) + sqrt(Gamma *pressure / density);
+	if (tmp>CFL*dx / dtx)
+	{
+		dtx = CFL*dx / tmp;
+	}
+	if (dtx > 0.25*CFL*dx*dx / Mu&&tau_type == NS&&Mu>0)
+	{
+		dtx = 0.25*CFL*dx*dx / Mu;
+	}
+	return dtx;
+}
+

code/function.h

@@ -0,0 +1,53 @@
+#pragma once
+#include<iostream>
+#include<cmath>
+#include<assert.h>
+using namespace std;
+
+extern int K;
+extern double Gamma;
+extern double Mu;
+extern double Nu;
+extern double T_inf; //the real temperature
+extern double R_gas; //the gas constant for real ideal air
+enum TAU_TYPE { Euler, NS, Sutherland, Power_law };
+extern TAU_TYPE tau_type;
+
+
+void Copy_Array(double* target, double* origin, int dim);
+void Array_zero(double* target, int dim);
+void Array_subtract(double* target, double* subtractor, int dim);
+void Array_scale(double* target, int dim, double value);
+
+void Convar_to_primvar_1D(double* primvar, double convar[3]);
+double Pressure(double density, double densityu, double densityE);
+double Temperature(double density, double pressure);
+double entropy(double density, double pressure);
+double Soundspeed(double density, double pressure);
+void Primvar_to_convar_1D(double* convar, double primvar[3]);
+double DensityU(double density, double u);
+double DensityE(double density, double u, double pressure);
+void Convar_to_char1D(double* character, double primvar[3], double convar[3]);
+void Char_to_convar1D(double* convar, double primvar[3], double charvar[3]);
+
+double U(double density, double q_densityu);
+double V(double density, double q_densityv);
+double Pressure(double density, double q_densityu, double q_densityv, double q_densityE);
+
+double Q_densityu(double density, double u);
+double Q_densityv(double density, double v);
+double Q_densityE(double density, double u, double v, double pressure);
+
+
+void Primvar_to_convar_2D(double* convar, double primvar[4]);//2d version
+void Convar_to_primvar_2D(double* primvar, double convar[4]);//2d version
+
+void Convar_to_char(double* character, double* primvar, double convar[4]);
+void Char_to_convar(double* convar, double* primvar, double character[4]);
+
+void Global_to_Local(double* change, double* normal);
+void Global_to_Local(double* change, double* origin, double* normal);
+void Local_to_Global(double* change, double* normal);
+void Local_to_Global(double* change, double* origin, double* normal);
+
+double Dtx(double dtx, double dx, double CFL, double density, double u, double v, double pressure);

code/fvm_gks1d.h

@@ -0,0 +1,158 @@
+#pragma once
+#include"function.h"
+#include<omp.h>
+#include"gks_basic.h"
+//extern type variables must be initialized (could not give value) in the corresponding .cpp file
+enum GKS1d_type{nothing, kfvs1st, kfvs2nd, gks1st, gks2nd};
+extern GKS1d_type gks1dsolver;
+enum Reconstruction_variable{conservative,characteristic};
+extern Reconstruction_variable reconstruction_variable;
+enum WENOtype { linear, wenojs, wenoz };
+extern WENOtype wenotype;
+extern bool is_reduce_order_warning;
+
+using namespace std;
+
+// to remember the geometry information
+typedef class Block1d
+{
+public:
+	//bool uniform;
+	int ghost;
+	int nodex; // mesh number
+	int nx;
+	int nodex_begin;
+	int nodex_end;
+	double dx;
+	double left;
+	double right;
+	int stages;
+	double timecoefficient[5][5][2];
+	double t; //current simulation time
+	double CFL; //cfl number, actually just a amplitude factor
+	double dt;//current_global time size
+	int step; //start from which step
+};
+
+// to remeber the cell avg values
+typedef class Fluid1d
+{
+public:
+	double primvar[3];
+	double convar[3];
+	double convar_old[3];
+	double cx; //center coordinate in x direction
+	double dx; //the mesh size dx
+};
+
+// to remember the fluid information in a fixed point,
+// such as reconstructioned value, or middle point value
+typedef class Point1d
+{
+public:
+	double convar[3];
+	double convar_old[3]; //this one is for compact point-wise reconstruction
+	double der1[3];
+	double x; // coordinate of a point
+};
+
+// remember the flux in a fixed interface,
+// for RK method, we only need f
+// for 2nd der method, we need derf
+typedef class Flux1d
+{
+public:
+	double F[3]; //total flux in dt time
+	double f[3]; // the f0 in t=0
+	double derf[4]; // the f_t in t=0
+};
+
+// every interfaces have left center and right value.
+// and center value for GKS
+// and a group of flux for RK method or multi-stage GKS method
+typedef class Interface1d
+{
+public:
+	Point1d left;
+	Point1d center;
+	Point1d right;
+	Flux1d* flux;
+	double x; // coordinate of the interface, equal to point1d.x
+};
+
+void Convar_to_primvar(Fluid1d *fluids, Block1d block);
+
+typedef void(*BoundaryCondition) (Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+void free_boundary_left(Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+void free_boundary_right(Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+void reflection_boundary_left(Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+void reflection_boundary_right(Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+void periodic_boundary_left(Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+void periodic_boundary_right(Fluid1d *fluids, Block1d block, Fluid1d bcvalue);
+
+
+void Reconstruction_within_cell(Interface1d *interfaces, Fluid1d *fluids, Block1d block);
+typedef void(*Reconstruction_within_Cell)(Point1d &left, Point1d &right, Fluid1d *fluids);
+extern Reconstruction_within_Cell cellreconstruction;
+
+void Check_Order_Reduce(Point1d &left, Point1d &right, Fluid1d &fluid);
+void Check_Order_Reduce_by_Lambda_1D(bool &order_reduce, double *convar);
+
+void Vanleer(Point1d &left, Point1d &right, Fluid1d *fluids);
+void weno_3rd_left(double& var, double& der1, double& der2,
+	double wn1, double w0, double wp1, double h);
+void weno_3rd_right(double& var, double& der1, double& der2,
+	double wn1, double w0, double wp1, double h);
+void WENO3_left(double& left, double w_back, double w, double w_forward, double h);
+void WENO3_right(double& right, double w_back, double w, double w_forward, double h);
+
+void WENO5_AO(Point1d &left, Point1d &right, Fluid1d *fluids);
+void weno_5th_ao_left(double &var, double &der1, double &der2, double wn2, double wn1, double w0, double wp1, double wp2, double h);
+void weno_5th_ao_right(double &var, double &der1, double &der2,	double wn2, double wn1, double w0, double wp1, double wp2,	double h);
+
+
+void Reconstruction_forg0(Interface1d *interfaces, Fluid1d *fluids, Block1d block);
+typedef void (*Reconstruction_forG0)(Interface1d &interfaces, Fluid1d *fluids);
+extern Reconstruction_forG0 g0reconstruction;
+void Center_2nd_collisionless(Interface1d& interfaces, Fluid1d* fluids);
+void Center_3rd(Interface1d& interfaces, Fluid1d* fluids);
+void Center_simple_avg(Interface1d &interfaces, Fluid1d *fluids);
+void Center_collision(Interface1d &interfaces, Fluid1d *fluids);
+
+
+void Calculate_flux(Flux1d** fluxes, Interface1d* interfaces, Block1d &block, int stage);
+typedef void(*Flux_function)(Flux1d &flux, Interface1d& interface, double dt);
+extern Flux_function flux_function;
+void GKS(Flux1d &flux, Interface1d& interface, double dt);
+
+
+typedef void(*TimeMarchingCoefficient)(Block1d &block);
+extern TimeMarchingCoefficient timecoe_list;
+void S1O1(Block1d& block);
+void S1O2(Block1d& block);
+void S2O4(Block1d& block);
+void RK4(Block1d &block);
+
+void Initial_stages(Block1d &block);
+Flux1d** Setflux_array(Block1d block);
+void SetUniformMesh(Block1d block, Fluid1d* fluids, Interface1d *interfaces, Flux1d **fluxes);
+
+
+void CopyFluid_new_to_old(Fluid1d *fluids, Block1d block);
+double Get_CFL(Block1d &block, Fluid1d *fluids, double tstop);
+double Dtx(double dtx, double dx, double CFL, double convar[3]);
+
+
+void Update(Fluid1d *fluids, Flux1d **fluxes, Block1d block, int stage);
+
+
+// two dimensional hllc and its subprogram
+void ESTIME(double& SL, double& SM, double& SR, double DL, double UL, double PL, double CL, double DR, double UR, double PR, double CR);
+void get_flux(double p[4], double* flux);
+void ustarforHLLC(double d1, double u1, double v1, double p1, double s1, double star1, double* ustar);
+void HLLC(Flux1d& flux, Interface1d& interface, double dt);
+void get_Euler_flux(double p[3], double* flux);
+void ustarforHLLC(double d1, double u1, double p1, double s1, double star1, double* ustar);
+
+void LF(Flux1d& flux, Interface1d& interface, double dt);
+

code/fvm_gks2d.h

@@ -0,0 +1,335 @@
+#pragma once
+#include"fvm_gks1d.h"
+extern int gausspoint;
+extern double *gauss_loc;
+extern double *gauss_weight;
+enum GKS2d_type { nothing_2d, kfvs1st_2d, kfvs2nd_2d, gks1st_2d, gks2nd_2d};
+extern GKS2d_type gks2dsolver;
+
+
+
+typedef class Block2d
+{
+	//geometry information of 2D block
+public:
+	int index;
+	bool uniform;
+	int ghost;
+	int nodex;
+	int nodey;
+	int nx; //total = real + ghost
+	int ny; //total = real + ghost
+	int xcell_begin;
+	int xcell_end;
+	int ycell_begin;
+	int ycell_end;
+	int xinterface_begin_n;
+	int xinterface_begin_t;
+	int xinterface_end_n;
+	int xinterface_end_t;
+	int yinterface_begin_n;
+	int yinterface_begin_t;
+	int yinterface_end_n;
+	int yinterface_end_t;
+	double dx;
+	double dy;
+	double overdx;
+	double overdy;
+	double left;
+	double right;
+	double down;
+	double up;
+	int stages;
+	double timecoefficient[5][5][2];
+	double t; //current simulation time
+	double CFL; //cfl number, actually just a amplitude factor
+	double dt;//current_global time size
+	int step; //start from which step
+	double T; //stop time
+	int outputperstep;
+};
+
+// to remember the fluid information in a fixed point,
+// such as reconstructioned value, or middle point value
+class Point2d
+{
+	// Note : store values relating to one point
+public:
+	double convar[4];
+	double der1x[4];
+	double der1y[4];
+	double x;
+	double y;
+	double normal[2];
+	bool is_reduce_order;
+};
+
+// to remeber the cell avg values
+//typedef class Fluid2d
+typedef class Fluid2d
+{
+	// Note : store the cell-averaged values
+	// do not have the global geometry/block information
+public:
+	double primvar[4];
+	double convar[4];
+	double convar_old[4];
+	double res[4];
+	double xindex;
+	double yindex;
+	double coordx;
+	double coordy;
+	double dx;
+	double dy;
+	double area;
+	bool is_hweno;
+	double node[8];
+	bool boundary;
+};
+
+
+// remember the flux in a fixed interface,
+// for RK method, we only need f
+// for 2nd der method, we need derf
+class Flux2d
+{
+public:
+	double x[4];
+	double f[4];
+	double derf[4];
+	double length;
+	double normal[2];
+};
+
+// every interfaces have left center and right value.
+// and center value for GKS
+// and a group of flux for RK method or multi-stage GKS method
+class Recon2d
+{
+	// Note : reconstruction can obtain left, cenrer, and right value
+public:
+	Point2d left;
+	Point2d center;
+	Point2d right;
+	double x;
+	double y;
+	double normal[2];
+};
+typedef class Interface2d
+{
+	// Note : for one interface, reconstruction can be for line-averagend or gauss point
+	// These reconstructed values can be used for the calculation of the flux at gauss points (stored in Flux2d_gauss) of the interface
+public:
+
+	Recon2d line;
+	Recon2d* gauss;
+	double x;
+	double y;
+	double normal[2];
+	double length;
+};
+
+class Flux2d_gauss
+{
+	// Note : Flux2d store the variables relating to flux
+	// Flux2d_gauss, includes several gauss points
+	// Final, one interface should have one Flux2d_gauss
+public:
+	Flux2d* gauss;
+};
+
+void SetGuassPoint();
+double Get_CFL(Block2d &block, Fluid2d *fluids, double tstop);
+
+
+
+void A_point(double *a, double* der, double* prim);
+void G_address(int no_u, int no_v, int no_xi, double *psi, double a[4], MMDF& m);
+void GL_address(int no_u, int no_v, int no_xi, double *psi, double a[4], MMDF& m);
+void GR_address(int no_u, int no_v, int no_xi, double *psi, double a[4], MMDF& m);
+
+
+
+typedef void(*BoundaryCondition2dnew) (Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void free_boundary_left(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void free_boundary_right(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void free_boundary_down(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void free_boundary_up(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void periodic_boundary_left(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void periodic_boundary_right(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void periodic_boundary_down(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void periodic_boundary_up(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void noslip_adiabatic_boundary_left(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void noslip_adiabatic_boundary(double* w_ghost, double* w_inner, double* normal, double* local_wall_velocity);
+void noslip_adiabatic_boundary_right(Fluid2d* fluids, Interface2d* face, Block2d block, Fluid2d bcvalue);
+void noslip_adiabatic_boundary_down(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void noslip_adiabatic_boundary_up(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void reflection_boundary_up(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void reflection_boundary_down(Fluid2d *fluids, Block2d block, Fluid2d bcvalue);
+void reflection_boundary_right(Fluid2d* fluids, Interface2d* face, Block2d block, Fluid2d bcvalue);
+void reflection_boundary_left(Fluid2d *fluids, Interface2d* face, Block2d block, Fluid2d bcvalue);
+void reflection_boundary(double* ghost, double* inner, double* normal);
+void inflow_boundary_left(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void inflow_boundary_right(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void inflow_boundary_up(Fluid2d* fluids, Block2d block, Fluid2d bcvalue);
+void ICfor_uniform_2d(Fluid2d* fluid, double* prim, Block2d block);
+
+bool negative_density_or_pressure(double* primvar);
+void Convar_to_Primvar(Fluid2d *fluids, Block2d block);
+void Convar_to_primvar(Fluid2d *fluid, Block2d block);
+void YchangetoX(double *fluidtmp, double *fluid);
+void XchangetoY(double *fluidtmp, double *fluid);
+void CopyFluid_new_to_old(Fluid2d *fluids, Block2d block);
+
+
+
+void Reconstruction_within_cell(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+typedef void(*Reconstruction_within_Cell_2D_normal)(Interface2d &left, Interface2d &right, Interface2d &down, Interface2d &up, Fluid2d *fluids, Block2d block);
+extern Reconstruction_within_Cell_2D_normal cellreconstruction_2D_normal;
+void First_order_normal(Interface2d& left, Interface2d& right, Interface2d& down, Interface2d& up, Fluid2d* fluids, Block2d block);
+void first_order(Point2d& left, Point2d& right, double* normal_l, double* normal_r, double* w);
+void Vanleer_normal(Interface2d &left, Interface2d &right, Interface2d &down, Interface2d &up, Fluid2d *fluids, Block2d block);
+void Vanleer_non_uniform(Point2d& left, Point2d& right, double* normal_l, double* normal_r, double* wn1, double* w, double* wp1, double h);
+void Vanleer(Point2d &left, Point2d &right, double *wn1, double *w, double *wp1, double h);
+void WENO3_normal(Interface2d& left, Interface2d& right, Interface2d& down, Interface2d& up, Fluid2d* fluids, Block2d block);
+void WENO3(Point2d& left, Point2d& right, double* wn1, double* w, double* wp1, double h);
+void WENO(Point2d& left, Point2d& right, double* wn2, double* wn1, double* w, double* wp1, double* wp2, double h);
+void WENO5_normal(Interface2d& left, Interface2d& right, Interface2d& down, Interface2d& up, Fluid2d* fluids, Block2d block);
+void WENO5_left(double& var, double& der1, double& der2, double wn2, double wn1, double w, double wp1, double wp2, double h);
+void WENO5_right(double& var, double& der1, double& der2, double wn2, double wn1, double w, double wp1, double wp2, double h);
+
+void WENO5_AO_normal(Interface2d &left, Interface2d &right, Interface2d &down, Interface2d &up, Fluid2d *fluids, Block2d block);
+void WENO5_AO(Point2d &left, Point2d &right, double * wn2, double * wn1, double * w, double * wp1, double * wp2, double h);
+
+
+
+typedef void(*Reconstruction_within_Cell_2D_tangent)(Interface2d *left, Interface2d *right, Interface2d *down, Interface2d *up, Fluid2d *fluids, Block2d block);
+extern Reconstruction_within_Cell_2D_tangent cellreconstruction_2D_tangent;
+void First_order_tangent(Interface2d* left, Interface2d* right, Interface2d* down, Interface2d* up, Fluid2d* fluids, Block2d block);
+void first_order_tangent(Point2d& gauss, Point2d& w0);
+void Vanleer_tangent(Interface2d *left, Interface2d *right, Interface2d *down, Interface2d *up, Fluid2d *fluids, Block2d block);
+void Vanleer(Point2d& gauss, Point2d &wn1, Point2d &w0, Point2d &wp1, double h);
+void Vanleer_tangential(double *coe, double wn1, double w0, double wp1, double h);
+
+void Polynomial_3rd(double* coefficent, double pn1, double w0, double pp1, double h);
+
+void Polynomial_5th(double* coefficent, double wn2, double wn1, double w0, double wp1, double wp2, double h);
+double Value_Polynomial(int order, double x0, double* coefficient);
+double Der1_Polynomial(int order, double x0, double* coefficient);
+
+void Polynomial_3rd_avg(double* coefficent, double wn1, double w0, double wp1, double h);
+double Value_Polynomial_3rd(double x, double coefficient[3]);
+double Der1_Polynomial_3rd(double x, double coefficient[3]);
+
+
+void WENO5_left(double& var, double wn2, double wn1, double w, double wp1, double wp2, double h);
+void WENO5_right(double& var, double wn2, double wn1, double w, double wp1, double wp2, double h);
+void WENO5_tangential(double* left, double* right, double* wn2, double* wn1, double* w, double* wp1, double* wp2, double h);
+void WENO5_tangential_for_slope(double* left, double* right, double* wn2, double* wn1, double* w, double* wp1, double* wp2, double h);
+
+void WENO5_tangential(Recon2d* re, Recon2d& wn2, Recon2d& wn1, Recon2d& w0, Recon2d& wp1, Recon2d& wp2, double h);
+void WENO5_tangent(Interface2d* left, Interface2d* right, Interface2d* down, Interface2d* up, Fluid2d* fluids, Block2d block);
+
+void WENO5_AO_tangent(Interface2d *left, Interface2d *right, Interface2d *down, Interface2d *up, Fluid2d *fluids, Block2d block);
+void WENO5_AO_tangential(Recon2d *re, Recon2d &wn2, Recon2d &wn1, Recon2d &w0, Recon2d &wp1, Recon2d &wp2, double h);
+void weno_5th_ao_2gauss(double &g1, double &g1x, double &g1xx, double &g2, double &g2x, double &g2xx, double wn2, double wn1, double w0, double wp1, double wp2, double h, int order);
+void WENO3_tangent(Interface2d* left, Interface2d* right, Interface2d* down, Interface2d* up, Fluid2d* fluids, Block2d block);
+void WENO3_tangential(Recon2d* re, Recon2d& wn1, Recon2d& w0, Recon2d& wp1, double h);
+void WENO3_tangential(double* left, double* right, double* wn1, double* w, double* wp1, double h);
+void WENO3_tangential_for_slope(double* left, double* right, double* wn1, double* w, double* wp1, double h);
+
+void Reconstruction_forg0(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+typedef void(*Reconstruction_forG0_2D_normal)(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+extern Reconstruction_forG0_2D_normal g0reconstruction_2D_normal;
+void Center_3rd_normal(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+void Center_3rd_Splitting(Interface2d& interface, double* w, double* wp1, double h);
+
+void Center_collision_normal(Interface2d* xinterfaces, Interface2d* yinterfaces, Fluid2d* fluids, Block2d block);
+void Center_do_nothing_normal(Interface2d* xinterfaces, Interface2d* yinterfaces, Fluid2d* fluids, Block2d block);
+
+void Center_4th_normal(Interface2d* xinterfaces, Interface2d* yinterfaces, Fluid2d* fluids, Block2d block);
+void Center_4th_normal(Interface2d& interface, double* wn1, double* w, double* wp1, double* wp2, double h);
+
+void Center_5th_normal(Interface2d* xinterfaces, Interface2d* yinterfaces, Fluid2d* fluids, Block2d block);
+void Center_5th_normal(Interface2d& interface, double* wn1, double* w, double* wp1, double* wp2, double h);
+
+
+
+typedef void(*Reconstruction_forG0_2D_tangent)(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+extern Reconstruction_forG0_2D_tangent g0reconstruction_2D_tangent;
+void Center_all_collision_multi(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+void Center_all_collision_2d_multi(Recon2d &gauss);
+void Center_first_order_tangent(Interface2d* xinterfaces, Interface2d* yinterfaces, Fluid2d* fluids, Block2d block);
+void Center_3rd_tangent(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+void Center_3rd_Multi(Point2d& gauss, Point2d& wn1, Point2d& w0, Point2d& wp1, double h);
+
+void Center_5th_tangent(Interface2d *xinterfaces, Interface2d *yinterfaces, Fluid2d *fluids, Block2d block);
+void Center_5th_Multi(Recon2d *re, Recon2d &wn2, Recon2d &wn1, Recon2d &w0, Recon2d &wp1, Recon2d &wp2, double h);
+
+
+
+
+void Calculate_flux(Flux2d_gauss** xfluxes, Flux2d_gauss** yfluxes, Interface2d* xinterfaces, Interface2d* yinterfaces, Block2d block, int stage);
+typedef void(*Flux_function_2d)(Flux2d &flux, Recon2d & re, double dt);
+extern Flux_function_2d flux_function_2d;
+void GKS2D_smooth(Flux2d &flux, Recon2d& interface, double dt);
+void GKS2D(Flux2d &flux, Recon2d& interface, double dt);
+void LF2D(Flux2d &flux, Recon2d& interface, double dt);
+void HLLC2D(Flux2d& flux, Recon2d& interface, double dt);
+void NS_by_central_difference_prim_2D(Flux2d& flux, Recon2d& interface, double dt);
+void NS_by_central_difference_convar_2D(Flux2d& flux, Recon2d& interface, double dt);
+
+
+
+typedef void(*TimeMarchingCoefficient_2d)(Block2d &block);
+extern TimeMarchingCoefficient_2d timecoe_list_2d;
+
+
+void S1O1_2D(Block2d &block);
+void S1O2_2D(Block2d& block);
+void RK2_2D (Block2d& block);
+void S2O4_2D(Block2d &block);
+void RK4_2D(Block2d& block);
+void Initial_stages(Block2d &block);
+
+
+
+Fluid2d *Setfluid(Block2d &block);
+
+Interface2d *Setinterface_array(Block2d block);
+
+Flux2d_gauss** Setflux_gauss_array(Block2d block);
+
+
+
+void SetUniformMesh(Block2d& block, Fluid2d* fluids, Interface2d *xinterfaces, Interface2d *yinterfaces, Flux2d_gauss **xfluxes, Flux2d_gauss **yfluxes);
+void Copy_geo_from_interface_to_line(Interface2d& interface);
+void Copy_geo_from_interface_to_flux(Interface2d& interface, Flux2d_gauss* flux, int stages);
+void Set_Gauss_Coordinate(Interface2d &xinterface, Interface2d &yinterface );
+void Set_Gauss_Intergation_Location_x(Point2d &xgauss, int index, double h);
+void Set_Gauss_Intergation_Location_y(Point2d &ygauss, int index, double h);
+
+
+
+
+
+void Update(Fluid2d *fluids, Flux2d_gauss **xfluxes, Flux2d_gauss **yfluxes, Block2d block, int stage);
+void Update_with_gauss(Fluid2d *fluids, Flux2d_gauss **xfluxes, Flux2d_gauss **yfluxes, Block2d block, int stage);
+
+void Check_Order_Reduce_by_Lambda_2D(bool &order_reduce, double *convar);
+void Recompute_KFVS_1st(Fluid2d &fluid,	double* center, double* left, double* right, double* down, double* up, Block2d& block);
+void KFVS_1st(double *flux, double* left, double* right, double dt);
+
+void SetNonUniformMesh(Block2d& block, Fluid2d* fluids,
+	Interface2d* xinterfaces, Interface2d* yinterfaces,
+	Flux2d_gauss** xfluxes, Flux2d_gauss** yfluxes, double** node2d);
+void Set_cell_geo_from_quad_node
+(Fluid2d& fluid, double* n1, double* n2, double* n3, double* n4);
+void Set_interface_geo_from_two_node
+(Interface2d& interface, double* n1, double* n2, int direction);
+void Set_Gauss_Coordinate_general_mesh_x
+(Interface2d& xinterface, double* node0, double* nodex1);
+void Set_Gauss_Coordinate_general_mesh_y
+(Interface2d& yinterface, double* node0, double* nodey1);
+
+void Residual2d(Fluid2d* fluids, Block2d block, int outputstep);

code/gks_basic.cpp

@@ -0,0 +1,520 @@
+//this libarary includes some basic functions for gks
+#include"gks_basic.h"
+#define S 110.4
+
+//artifical viscosity for euler problem
+double c1_euler = -1.0; //initialization
+double c2_euler = -1.0; //initialization
+bool is_Prandtl_fix = false; //initialization
+double Pr = 1.0; //initialization
+
+
+
+// to prepare the basic element for moment calculation
+MMDF1d::MMDF1d() { u = 1.0; lambda = 1.0; };
+MMDF1d::MMDF1d(double u_in, double lambda_in)
+{
+	u = u_in;
+	lambda = lambda_in;
+	calcualte_MMDF1d();
+}
+void MMDF1d::calcualte_MMDF1d()
+{
+	uwhole[0] = 1;
+	uwhole[1] = u;
+	uplus[0] = 0.5 * Alpha(lambda, -u);
+	uminus[0] = 0.5 * Alpha(lambda, u);
+	uplus[1] = u * uplus[0] + 0.5 * Beta(lambda, u);
+	uminus[1] = u * uminus[0] - 0.5 * Beta(lambda, u);
+	for (int i = 2; i <= 9; i++)
+	{
+		uwhole[i] = u * uwhole[i - 1] + 0.5 * (i - 1) / lambda * uwhole[i - 2];
+	}
+	for (int i = 2; i <= 9; i++)
+	{
+		uplus[i] = u * uplus[i - 1] + 0.5 * (i - 1) / lambda * uplus[i - 2];
+		uminus[i] = u * uminus[i - 1] + 0.5 * (i - 1) / lambda * uminus[i - 2];
+	}
+	xi2 = 0.5 * K / lambda;
+	xi4 = 0.25 * (K * K + 2 * K) / (lambda * lambda);
+	//xi6?? how to calculate
+	xi6 = 0.5 * (K + 4) / lambda * xi4;
+
+	for (int i = 0; i < 10; i++)
+	{
+		for (int k = 0; k < 4; k++)
+		{
+			if ((i + 2 * k) <= 9)
+			{
+				if (k == 0)
+				{
+					upxi[i][k] = uplus[i];
+					unxi[i][k] = uminus[i];
+				}
+				if (k == 1)
+				{
+					upxi[i][k] = uplus[i] * xi2;
+					unxi[i][k] = uminus[i] * xi2;
+				}
+				if (k == 2)
+				{
+					upxi[i][k] = uplus[i] * xi4;
+					unxi[i][k] = uminus[i] * xi4;
+				}
+				if (k == 3)
+				{
+					upxi[i][k] = uplus[i] * xi6;
+					unxi[i][k] = uminus[i] * xi6;
+				}
+			}
+		}
+	}
+	for (int i = 0; i < 10; i++)
+	{
+		for (int k = 0; k < 4; k++)
+		{
+			if ((i + 2 * k) <= 9)
+			{
+				if (k == 0)
+				{
+					uxi[i][k] = uwhole[i];
+				}
+				if (k == 1)
+				{
+					uxi[i][k] = uwhole[i] * xi2;
+				}
+				if (k == 2)
+				{
+					uxi[i][k] = uwhole[i] * xi4;
+				}
+				if (k == 3)
+				{
+					uxi[i][k] = uwhole[i] * xi6;
+				}
+			}
+		}
+	}
+}
+
+//a general G function
+void G(int no_u, int no_xi, double* psi, double a[3], MMDF1d m)
+{
+	psi[0] = a[0] * m.uxi[no_u][no_xi] + a[1] * m.uxi[no_u + 1][no_xi] + a[2] * 0.5 * (m.uxi[no_u + 2][no_xi] + m.uxi[no_u][no_xi + 1]);
+	psi[1] = a[0] * m.uxi[no_u + 1][no_xi] + a[1] * m.uxi[no_u + 2][no_xi] + a[2] * 0.5 * (m.uxi[no_u + 3][no_xi] + m.uxi[no_u + 1][no_xi + 1]);
+	psi[2] = 0.5 * (a[0] * (m.uxi[no_u + 2][no_xi] + m.uxi[no_u][no_xi + 1]) +
+		a[1] * (m.uxi[no_u + 3][no_xi] + m.uxi[no_u + 1][no_xi + 1]) +
+		a[2] * 0.5 * (m.uxi[no_u + 4][no_xi] + m.uxi[no_u][no_xi + 2] + 2 * m.uxi[no_u + 2][no_xi + 1]));
+
+}
+void GL(int no_u, int no_xi, double* psi, double a[3], MMDF1d m)
+{
+	psi[0] = a[0] * m.upxi[no_u][no_xi] + a[1] * m.upxi[no_u + 1][no_xi] + a[2] * 0.5 * (m.upxi[no_u + 2][no_xi] + m.upxi[no_u][no_xi + 1]);
+	psi[1] = a[0] * m.upxi[no_u + 1][no_xi] + a[1] * m.upxi[no_u + 2][no_xi] + a[2] * 0.5 * (m.upxi[no_u + 3][no_xi] + m.upxi[no_u + 1][no_xi + 1]);
+	psi[2] = 0.5 * (a[0] * (m.upxi[no_u + 2][no_xi] + m.upxi[no_u][no_xi + 1]) +
+		a[1] * (m.upxi[no_u + 3][no_xi] + m.upxi[no_u + 1][no_xi + 1]) +
+		a[2] * 0.5 * (m.upxi[no_u + 4][no_xi] + m.upxi[no_u][no_xi + 2] + 2 * m.upxi[no_u + 2][no_xi + 1]));
+
+}
+void GR(int no_u, int no_xi, double* psi, double a[3], MMDF1d m)
+{
+	psi[0] = a[0] * m.unxi[no_u][no_xi] + a[1] * m.unxi[no_u + 1][no_xi] + a[2] * 0.5 * (m.unxi[no_u + 2][no_xi] + m.unxi[no_u][no_xi + 1]);
+	psi[1] = a[0] * m.unxi[no_u + 1][no_xi] + a[1] * m.unxi[no_u + 2][no_xi] + a[2] * 0.5 * (m.unxi[no_u + 3][no_xi] + m.unxi[no_u + 1][no_xi + 1]);
+	psi[2] = 0.5 * (a[0] * (m.unxi[no_u + 2][no_xi] + m.unxi[no_u][no_xi + 1]) +
+		a[1] * (m.unxi[no_u + 3][no_xi] + m.unxi[no_u + 1][no_xi + 1]) +
+		a[2] * 0.5 * (m.unxi[no_u + 4][no_xi] + m.unxi[no_u][no_xi + 2] + 2 * m.unxi[no_u + 2][no_xi + 1]));
+
+}
+
+//moments of the maxwellian distribution function
+MMDF::MMDF() { u = 1.0; v = 1.0; lambda = 1.0; };
+MMDF::MMDF(double u_in, double v_in, double lambda_in)
+{
+	u = u_in;
+	v = v_in;
+	lambda = lambda_in;
+	calcualte_MMDF();
+}
+void MMDF::calcualte_MMDF()
+{
+	//the moment of u in whole domain
+	uwhole[0] = 1;
+	uwhole[1] = u;
+	//the moment of u in half domain (>0 && <0)
+	uplus[0] = 0.5 * Alpha(lambda, -u);
+	uminus[0] = 1.0 - uplus[0];
+	uplus[1] = u * uplus[0] + 0.5 * Beta(lambda, u);
+	uminus[1] = u - uplus[1];
+	double overlambda = 1.0 / lambda;
+	for (int i = 2; i <= 6; i++)
+	{
+		uwhole[i] = u * uwhole[i - 1] + 0.5 * (i - 1) * overlambda * uwhole[i - 2];
+	}
+	for (int i = 2; i <= 6; i++)
+	{
+		uplus[i] = u * uplus[i - 1] + 0.5 * (i - 1) * overlambda * uplus[i - 2];
+		uminus[i] = uwhole[i] - uplus[i];
+	}
+	//the moment of v in whole domain
+	//the moment of v in half domain is useless
+	vwhole[0] = 1;
+	vwhole[1] = v;
+	for (int i = 2; i <= 6; i++)
+	{
+		vwhole[i] = v * vwhole[i - 1] + 0.5 * (i - 1) * overlambda * vwhole[i - 2];
+	}
+	//the moment of xi (kesi)
+	xi2 = 0.5 * K * overlambda;
+	xi4 = 0.25 * (K * K + 2 * K) * overlambda * overlambda;
+	//xi6?? how to calculate
+	//xi6 = 0.5*(K + 4) / lambda*xi4; // no use
+
+	//get other variables by the above mement variables
+	//i, j, k, means i power of u, j power of v, and 2k power of xi
+	for (int i = 0; i < 7; i++)
+	{
+		for (int j = 0; j < 7; j++)
+		{
+			for (int k = 0; k < 3; k++)
+			{
+				if ((i + j + 2 * k) <= 6)
+				{
+					if (k == 0)
+					{
+						// u, half domain > 0; v, whole domain; xi, equals 0
+						upvxi[i][j][k] = uplus[i] * vwhole[j];
+						// u, half domain < 0; v, whole domain; xi, equals 0
+						unvxi[i][j][k] = uminus[i] * vwhole[j];
+						// u, whole domain ; v, whole domain; xi, equals 0
+						uvxi[i][j][k] = upvxi[i][j][k] + unvxi[i][j][k];
+						// = (uplus[i] + uminus[i]) * vwhole[j] = uwhole[i] * vwhole[j]
+					}
+					if (k == 1)
+					{
+						// u, half domain > 0; v, whole domain; xi, equals 2
+						upvxi[i][j][k] = uplus[i] * vwhole[j] * xi2;
+						// u, half domain < 0; v, whole domain; xi, equals 2
+						unvxi[i][j][k] = uminus[i] * vwhole[j] * xi2;
+						// u, whole domain ; v, whole domain; xi, equals 2
+						uvxi[i][j][k] = upvxi[i][j][k] + unvxi[i][j][k];
+					}
+					if (k == 2)
+					{
+						// u, half domain > 0; v, whole domain; xi, equals 4
+						upvxi[i][j][k] = uplus[i] * vwhole[j] * xi4;
+						// u, half domain < 0; v, whole domain; xi, equals 4
+						unvxi[i][j][k] = uminus[i] * vwhole[j] * xi4;
+						// u, whole domain ; v, whole domain; xi, equals 4
+						uvxi[i][j][k] = upvxi[i][j][k] + unvxi[i][j][k];
+					}
+
+				}
+			}
+		}
+	}
+}
+
+MMDF1st::MMDF1st() { u = 1.0; v = 1.0; lambda = 1.0; };
+MMDF1st::MMDF1st(double u_in, double v_in, double lambda_in)
+{
+	u = u_in;
+	v = v_in;
+	lambda = lambda_in;
+	calcualte_MMDF1st();
+}
+void MMDF1st::calcualte_MMDF1st()
+{
+	uwhole[0] = 1.0;
+	uwhole[1] = u;
+
+	//return erfc(sqrt(lambda)*u);
+
+	//return exp(-lambda*u*u) / sqrt(pi*lambda);
+	uminus[0] = 0.5 * Alpha(lambda, u);
+	uplus[0] = 1.0 - uminus[0];
+	double beta = Beta(lambda, u);
+	//uplus[1] = u*uplus[0] + 0.5*Beta(lambda, u);
+	//uminus[1] = u*uminus[0] - 0.5*Beta(lambda, u);
+	uplus[1] = u * uplus[0] + 0.5 * beta;
+	uminus[1] = u * uminus[0] - 0.5 * beta;
+	double overlambda = 1.0 / lambda;
+	for (int i = 2; i <= 3; i++)
+	{
+		uwhole[i] = u * uwhole[i - 1] + 0.5 * (i - 1) * overlambda * uwhole[i - 2];
+	}
+	for (int i = 2; i <= 3; i++)
+	{
+		uplus[i] = u * uplus[i - 1] + 0.5 * (i - 1) * overlambda * uplus[i - 2];
+		uminus[i] = uwhole[i] - uplus[i];
+	}
+	vwhole[0] = 1;
+	vwhole[1] = v;
+	for (int i = 2; i <= 2; i++)
+	{
+		vwhole[i] = v * vwhole[i - 1] + 0.5 * (i - 1) * overlambda * vwhole[i - 2];
+	}
+	xi2 = 0.5 * K * overlambda;
+}
+
+
+
+void Collision(double* w0, double left, double right, MMDF& m2, MMDF& m3)
+{
+	// get the equilibrium variables by collision
+	w0[0] = left * m2.uplus[0] + right * m3.uminus[0];
+	w0[1] = left * m2.uplus[1] + right * m3.uminus[1];
+	w0[2] = left * m2.vwhole[1] * m2.uplus[0] + right * m3.vwhole[1] * m3.uminus[0];
+	w0[3] = 0.5 * left * (m2.uplus[2] + m2.uplus[0] * m2.vwhole[2] + m2.uplus[0] * m2.xi2) +
+		0.5 * right * (m3.uminus[2] + m3.uminus[0] * m3.vwhole[2] + m3.uminus[0] * m3.xi2);
+}
+void Collision(double* w0, double left, double right, MMDF1st& m2, MMDF1st& m3)
+{
+	w0[0] = left * m2.uplus[0] + right * m3.uminus[0];
+	w0[1] = left * m2.uplus[1] + right * m3.uminus[1];
+	w0[2] = left * m2.vwhole[1] * m2.uplus[0] + right * m3.vwhole[1] * m3.uminus[0];
+	w0[3] = 0.5 * left * (m2.uplus[2] + m2.uplus[0] * m2.vwhole[2] + m2.uplus[0] * m2.xi2) +
+		0.5 * right * (m3.uminus[2] + m3.uminus[0] * m3.vwhole[2] + m3.uminus[0] * m3.xi2);
+}
+
+
+
+void A(double* a, double der[4], double prim[4])
+{
+	//Note: solve the a coefficient (slope information), equivalent to solving the matrix to get a
+	//input: der[4], the gradient of conservative variables, density, momentumX, momentumY, energy, partialW/patialX
+	//input: prim[4], the primary variables (density, u, v, lambda) at one point or one interface
+	//output: a, solve the results
+
+	double R4, R3, R2;
+	double overden = 1.0 / prim[0];
+	//double overlambda = 1.0 / prim[3];
+	R4 = der[3] * overden - 0.5 * (prim[1] * prim[1] + prim[2] * prim[2] + 0.5 * (K + 2) / prim[3]) * der[0] * overden;
+	R3 = (der[2] - prim[2] * der[0]) * overden;
+	R2 = (der[1] - prim[1] * der[0]) * overden;
+	a[3] = (4.0 / (K + 2)) * prim[3] * prim[3] * (2 * R4 - 2 * prim[1] * R2 - 2 * prim[2] * R3);
+	a[2] = 2 * prim[3] * R3 - prim[2] * a[3];
+	a[1] = 2 * prim[3] * R2 - prim[1] * a[3];
+	a[0] = der[0] * overden - prim[1] * a[1] - prim[2] * a[2] - 0.5 * a[3] * (prim[1] * prim[1] + prim[2] * prim[2] + 0.5 * (K + 2) / prim[3]);
+
+	//R4 = der[3] / prim[0] - 0.5*(prim[1] * prim[1] + prim[2] * prim[2] + 0.5*(K + 2) / prim[3])*der[0] / prim[0];
+	//R3 = (der[2] - prim[2]*der[0]) / prim[0];
+	//R2 = (der[1] - prim[1]*der[0]) / prim[0];
+	//a[3] = 4.0 * prim[3]*prim[3] / (K + 2)*(2.0 * R4 - 2.0 * prim[1]*R2 - 2.0 * prim[2]*R3);
+	//a[2] = 2.0 * prim[3]*R3 - prim[2]*a[3];
+	//a[1] = 2.0 * prim[3]*R2 - prim[1]*a[3];
+	//a[0] = der[0] / prim[0] - prim[1]*a[1] - prim[2]*a[2] - 0.5*a[3] * (prim[1]*prim[1] + prim[2]*prim[2] + 0.5*(K + 2) / prim[3]);
+
+}
+
+
+
+double Get_Tau_NS(double density0, double lambda0)
+{
+	if (tau_type == Euler)
+	{
+		return 0.0;
+	}
+	else if (tau_type == NS)
+	{
+		if (Mu > 0.0)
+		{
+			//cout << "here" << endl;
+			return 2.0 * Mu * lambda0 / density0;
+		}
+		else if (Nu > 0.0)
+		{
+			return 2 * Nu * lambda0;
+		}
+		else
+		{
+			return 0.0;
+		}
+	}
+	else if (tau_type == Sutherland)
+	{
+		return TauNS_Sutherland(density0, lambda0);
+	}
+	else
+	{
+		return TauNS_power_law(density0, lambda0);
+	}
+
+}
+double Get_Tau(double density_left, double density_right, double density0, double lambda_left, double lambda_right, double lambda0, double dt)
+{
+	if (tau_type == Euler)
+	{
+		if (c1_euler <= 0 && c2_euler <= 0)
+		{
+			return 0.0;
+		}
+		else
+		{
+			double C = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right);
+			//C+= c2_euler * abs(density_left - density_right) / abs(density_left + density_right);
+			//if (C < 10)
+			//{
+			return c1_euler * dt + dt * C;
+			//}
+			//else
+			//{
+			//	return c1_euler*dt + 10*dt;
+			//	//return c1_euler*dt + dt*C;
+			//}
+		}
+	}
+	else if (tau_type == NS)
+	{
+		double tau_n = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right) * dt;
+		if (tau_n != tau_n)
+		{
+			tau_n = 0.0;
+		}
+		if ((Mu > 0.0 && Nu > 0.0) || (Mu < 0.0 && Nu < 0.0))
+		{
+			return 0.0;
+		}
+		else
+		{
+			if (Mu > 0.0)
+			{
+				return tau_n + 2.0 * Mu * lambda0 / density0;
+			}
+			else if (Nu > 0.0)
+			{
+				return tau_n + 2 * Nu * lambda0;
+			}
+			else
+			{
+				return 0.0;
+			}
+		}
+	}
+	else if (tau_type == Sutherland)
+	{
+		double tau_n = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right) * dt;
+		if (tau_n != tau_n)
+		{
+			tau_n = 0.0;
+		}
+		if ((Mu > 0.0 && Nu > 0.0) || (Mu < 0.0 && Nu < 0.0))
+		{
+
+			return 0.0;
+		}
+		else
+		{
+			if (Mu > 0.0)
+			{
+				return tau_n + TauNS_Sutherland(density0, lambda0);
+
+			}
+			else if (Nu > 0.0)
+			{
+				return tau_n + 2 * Nu * lambda0;
+			}
+			else
+			{
+				return 0.0;
+			}
+		}
+	}
+	else if (tau_type == Power_law)
+	{
+
+		double tau_n = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right) * dt;
+		if (tau_n != tau_n)
+		{
+			tau_n = 0.0;
+		}
+		if ((Mu > 0.0 && Nu > 0.0) || (Mu < 0.0 && Nu < 0.0))
+		{
+			return 0.0;
+		}
+		else
+		{
+			if (Mu > 0.0)
+			{
+				return tau_n + TauNS_power_law(density0, lambda0);
+			}
+			else if (Nu > 0.0)
+			{
+				return tau_n + 2 * Nu * lambda0;
+			}
+			else
+			{
+				return 0.0;
+			}
+		}
+
+	}
+	else
+	{
+		return 0.0;
+	}
+}
+double TauNS_Sutherland(double density0, double lambda0)
+{
+	double T0 = 1.0 / 2.0 / (lambda0 * R_gas);
+	double mu = Mu * pow(T0 / T_inf, 1.5) * (T_inf + S) / (T0 + S);
+	return 2.0 * mu * lambda0 / density0;
+}
+double TauNS_power_law(double density0, double lambda0)
+{
+	double T0 = 1.0 / 2.0 / (lambda0 * R_gas);
+	double mu = Mu * pow(T0 / T_inf, 0.76);
+	return 2.0 * mu * lambda0 / density0;
+}
+
+
+
+//计算λ
+double Lambda(double density, double u, double densityE)
+{
+	return (K + 1.0) * 0.25 * (density / (densityE - 0.5 * density * (u * u)));
+}
+double Lambda(double density, double u, double v, double densityE)
+{
+	return (K + 2.0) * 0.25 * (density / (densityE - 0.5 * density * (u * u + v * v)));
+}
+
+void Convar_to_ULambda_1d(double* primvar, double convar[3])
+{
+	primvar[0] = convar[0];
+	primvar[1] = U(convar[0], convar[1]);
+	primvar[2] = Lambda(convar[0], primvar[1], convar[2]);
+}
+
+//计算α，α=erfc（±√λ*U）
+double Alpha(double lambda, double u)
+{
+	return erfc(sqrt(lambda) * u);
+	//return erfcmy(sqrt(lambda)*u);
+}
+//计算β，β=e^(-λ*U^2)/(√πλ)
+double Beta(double lambda, double u)
+{
+	double pi = 3.14159265358979323846;
+	return exp(-lambda * u * u) / sqrt(pi * lambda);
+}
+
+
+//solution of matrix equation b=Ma 1D
+void Microslope(double* a, double der[3], double prim[3])
+{
+	double R4, R2;
+	R4 = der[2] / prim[0] - 0.5 * (prim[1] * prim[1] + 0.5 * (K + 1) / prim[2]) * der[0] / prim[0];
+	R2 = (der[1] - prim[1] * der[0]) / prim[0];
+	a[2] = 4 * prim[2] * prim[2] / (K + 1) * (2 * R4 - 2 * prim[1] * R2);
+	a[1] = 2 * prim[2] * R2 - prim[1] * a[2];
+	a[0] = der[0] / prim[0] - prim[1] * a[1] - 0.5 * a[2] * (prim[1] * prim[1] + 0.5 * (K + 1) / prim[2]);
+}
+
+void Convar_to_ULambda_2d(double* primvar, double convar[4]) //density, U, V, lambda
+{
+	primvar[0] = convar[0];
+	primvar[1] = U(convar[0], convar[1]);
+	primvar[2] = V(convar[0], convar[2]);
+	//here primvar[3] refers lambda
+	primvar[3] = Lambda(convar[0], primvar[1], primvar[2], convar[3]);
+}
+

code/gks_basic.h

@@ -0,0 +1,105 @@
+#pragma once
+#include"function.h"
+
+
+extern double c1_euler;
+extern double c2_euler;
+extern bool is_Prandtl_fix;
+extern double Pr;
+
+
+
+//basic gks function
+// to store the moment
+class MMDF1d
+{
+private:
+	double u;
+	double lambda;
+
+public:
+	double uwhole[10];
+	double uplus[10];
+	double uminus[10];
+	double upxi[10][4];
+	double unxi[10][4];
+	double uxi[10][4];
+	double xi2;
+	double xi4;
+	double xi6;
+	MMDF1d();
+	MMDF1d(double u_in, double lambda_in);
+	void calcualte_MMDF1d();
+};
+
+// to calculate the microsolpe moment
+void G(int no_u, int no_xi, double* psi, double a[3], MMDF1d m);
+void GL(int no_u, int no_xi, double* psi, double a[3], MMDF1d m);
+void GR(int no_u, int no_xi, double* psi, double a[3], MMDF1d m);
+
+//moments of the maxwellian distribution function
+class MMDF
+{
+private:
+	double u;
+	double v;
+	double lambda;
+
+public:
+	double uwhole[7];
+	double uplus[7];
+	double uminus[7];
+	double vwhole[7];
+	double upvxi[7][7][3];
+	double unvxi[7][7][3];
+	double uvxi[7][7][3];
+	double xi2;
+	double xi4;
+	MMDF();
+	MMDF(double u_in, double v_in, double lambda_in);
+	void calcualte_MMDF();
+};
+
+class MMDF1st
+{
+private:
+	double u;
+	double v;
+	double lambda;
+
+public:
+	double uwhole[4];
+	double uplus[4];
+	double uminus[4];
+	double vwhole[3];
+	double xi2;
+	MMDF1st();
+	MMDF1st(double u_in, double v_in, double lambda_in);
+	void calcualte_MMDF1st();
+};
+
+
+void Collision(double *w0, double left, double right, MMDF &m2, MMDF &m3);
+void Collision(double *w0, double left, double right, MMDF1st &m2, MMDF1st &m3);
+
+void A(double *a, double der[4], double prim[4]);
+
+double Get_Tau_NS(double density0, double lambda0); // solve the smooth tau
+double TauNS_Sutherland(double density0, double lambda0); // solve the smooth tau by using sutherland
+double TauNS_power_law(double density0, double lambda0); //solver the smooth tau by using power law
+double Get_Tau(double density_left, double density_right, double density0, double lambda_left, double lambda_right, double lambda0, double dt); // solve not smooth numerical tau
+
+
+double Lambda(double density, double u, double densityE);
+double Lambda(double density, double u, double v, double densityE);
+void Convar_to_ULambda_1d(double* primvar, double convar[3]);
+
+//计算α，α=erfc（±√λ*U）
+double Alpha(double lambda, double u);
+//计算β，β=e^(-λ*U^2)/(√πλ)
+double Beta(double lambda, double u);
+//solution of matrix equation b=Ma
+void Microslope(double *a, double der[3], double prim[3]); //in one dimensional
+
+void Convar_to_ULambda_2d(double* primvar, double convar[4]);
+

code/main.cpp

@@ -0,0 +1,52 @@
+
+//说明指南:
+//这是一个气体动理学格式的求解器，支持二维单块结构化网格模拟，目的是方便有需求的科研工作者入门之用。
+//程序基本采用面向过程编写，没有什么弯弯绕，
+//故建议从main函数 顺序 看起，搭配注释学习。注释采用土洋结合，以简洁方便的传递作者意图为原则。
+//
+//如用作学术用途，请引用任意下列文献：
+//（1）X.JI, F.ZHAO, W.SHYY, & K.XU(2018).
+//A family of high - order gas - kinetic schemes and its comparison with Riemann solver based high - order methods.
+//Journal of Computational Physics, 356, 150 - 173.
+//（2）X.JI, & K.XU(2020).Performance Enhancement for High - order Gas - kinetic Scheme Based on WENO - adaptive - order Reconstruction.
+//Communication in Computational Physics, 28, 2, 539 - 590
+
+
+#include<omp.h>  //omp并行，尖括号<>代表外部库文件
+#include"output.h" //双引号代表该项目中的头文件
+
+//以下是几个测试算例的头文件
+#include"accuracy_test.h"
+#include"boundary_layer.h"
+#include"cylinder.h"
+#include"riemann_problem.h"
+//end
+
+using namespace std; //默认std的命名规则
+
+//程序运行开始，手动输入omp并行的线程数
+//如果在函数名上面注释，那么函数在别处调用时，可以看到注释
+void Set_omp_thread() 
+{
+	//在小于等于cpu核数的时候，n线程可以理解为n核并行
+	int num_thread;
+	cout << "please_specify threads number for omp parallel:  ";
+	cin >> num_thread;
+	omp_set_num_threads(num_thread);
+}
+
+int main()
+{
+	Set_omp_thread(); //设置omp并行线程
+	make_directory_for_result(); //兼容linux的结果文件夹的临时方案
+	//accuracy_sinwave_1d(); //一维精度测试，通过周期线性的正弦波传播算例，测试格式在光滑（smooth flow）无粘工况下的精度
+	//riemann_problem_1d();  //一维黎曼问题测试，测试格式对于可压缩间断问题的鲁棒性和分辨率
+    //accuracy_sinwave_2d(); //二维精度测试，通过二维周期线性的正弦波传播算例，测试格式在光滑（smooth flow）无粘工况下的精度
+	//riemann_problem_2d();  //二维黎曼问题测试，测试格式对于多维可压缩间断问题的鲁棒性和分辨率
+	boundary_layer();      //二维边界层问题测试，测试格式对于光滑（smooth flow）粘性边界的分辨率，非均匀直角网格
+	cylinder();              //二维超音速圆柱绕流问题测试，测试格式对于间断问题和非均匀网格的鲁棒性
+	
+    return 0;
+}
+
+

code/mesh_read.cpp

@@ -0,0 +1,331 @@
+#include"mesh_read.h"
+
+//兼容linux
+string add_mesh_directory_modify_for_linux()
+{
+#if defined(_WIN32)
+	return "";
+#else
+	return "../";
+#endif
+}
+
+void Read_structured_mesh(string& filename,
+	Fluid2d** fluids, Interface2d** xinterfaces, Interface2d** yinterfaces,
+	Flux2d_gauss*** xfluxes, Flux2d_gauss*** yfluxes, Block2d& block)
+{
+	//参考 因特网无名前辈
+	int block_num = 1;
+
+	ifstream in;
+	in.open(filename, ifstream::binary);
+	if (!in.is_open())
+	{
+		cout << "cannot find input structured grid file (.plt) " << endl;
+		while (1);
+		exit(0);
+	}
+	cout << "reading mesh..." << endl;
+
+	//prepare node information
+	double*** node2d = new double** [block_num];
+
+
+	// Header section
+	// section i
+	char buffer[9];
+	in.read(buffer, 8);
+	buffer[8] = '\0';
+	string version = string(buffer);
+	if (version != "#!TDV112")
+		cout << "Parser based on version #!TDV112, current version is: " << version << endl;
+
+	// section ii
+	int IntegerValue;
+	in.read((char*)&IntegerValue, sizeof(IntegerValue));
+	if (IntegerValue != 1)
+		cout << "the interger value is not equal to 1 but " << IntegerValue << endl;
+
+	// section iii
+	//file type  0=FULL 1=GRID 2=SOLUTION
+	int file_type;
+	in.read((char*)&file_type, sizeof(file_type));
+	//file title ---> defined name of the plt
+	string title;
+	read_ascii_to_string(title, in);
+	cout << "the mesh title is " << title;
+	//the number of total variables (x, y, z) as defult
+	int number_variables;
+	in.read((char*)&number_variables, sizeof(number_variables));
+	if (number_variables != 3)
+	{
+		cout << "warning, the input grid file has the number of variables " << number_variables << endl;
+	}
+	string variable_names[3];
+	for (int i = 0; i < number_variables; i++)
+		read_ascii_to_string(variable_names[i], in);
+
+	// section iv
+	// we try to read multi-zone; the marker for a start of new zone is a float number 
+	// validation_marker==299.0
+	float validation_marker;
+	in.read((char*)&validation_marker, sizeof(validation_marker));
+	assert(validation_marker == 299.0);
+	int zone_num = 0;
+	int imax[100];
+	int jmax[100];
+	int kmax[100];
+	size_t number_nodes[100];
+	while (validation_marker == 299.0)
+	{
+		zone_num++;
+		string zone_name;
+		int parent_zone;
+		int strand_id;
+		double solution_time;
+		int not_used;
+		int zone_type;  // 0=ORDERED 1=FELINESEG 2=FETRIANGLE 3=FEQUADRILATERAL 
+						// 4=FETETRAHEDRON 5=FEBRICK 6=FEPOLYGON 7=FEPOLYHEDRON
+		int data_packing; // 0=Block, 1=Point
+		int var_location;
+
+		read_ascii_to_string(zone_name, in);
+		in.read((char*)&parent_zone, sizeof(parent_zone));
+		in.read((char*)&strand_id, sizeof(strand_id));
+		in.read((char*)&solution_time, sizeof(solution_time));
+		in.read((char*)&not_used, sizeof(not_used));
+		in.read((char*)&zone_type, sizeof(zone_type));
+		in.read((char*)&data_packing, sizeof(data_packing));
+		in.read((char*)&var_location, sizeof(var_location));
+
+
+		if (var_location == 0)
+		{
+		}
+		else
+		{
+			cout << "we donnt support this function" << endl;
+			exit(0);
+		}
+
+		// face neighbors
+		int face_neighbors;
+		in.read((char*)&face_neighbors, sizeof(face_neighbors));
+
+		if (face_neighbors == 0)
+		{
+		}
+		else
+		{
+			cout << "we donnt support this function" << endl;
+			exit(0);
+		}
+
+		// ordered zone
+		if (zone_type == 0)
+		{
+			in.read((char*)&imax[zone_num - 1], sizeof(imax[zone_num - 1]));
+			in.read((char*)&jmax[zone_num - 1], sizeof(jmax[zone_num - 1]));
+			in.read((char*)&kmax[zone_num - 1], sizeof(kmax[zone_num - 1]));
+			number_nodes[zone_num - 1]
+				= imax[zone_num - 1] * jmax[zone_num - 1] * kmax[zone_num - 1];
+			cout << zone_name << " has total "
+				<< imax[zone_num - 1] << " * " << jmax[zone_num - 1] << " * " << kmax[zone_num - 1] << "points " << endl;
+
+		}
+		else
+		{
+			cout << "we donnt support unstructure mesh currently" << endl;
+			exit(0);
+		}
+		// to do create an auxiliary structure
+		// add all the structures to a vector
+		in.read((char*)&validation_marker, sizeof(validation_marker));
+
+
+		// read the next marker to check if 
+		// there is more than one zone
+		in.read((char*)&validation_marker, sizeof(validation_marker));
+	} // end of zone section
+
+	int sum_mesh_number = 0;
+
+	for (int iblock = 0; iblock < block_num; iblock++)
+	{
+		block.nodex = imax[iblock] - 1;
+		block.nodey = jmax[iblock] - 1;
+		for (int iblock = 0; iblock < block_num; iblock++)
+		{
+			block.nx = block.nodex + 2 * block.ghost;
+			block.ny = block.nodey + 2 * block.ghost;
+			sum_mesh_number += block.nodex * block.nodey;
+		}
+		node2d[iblock] = new double* [(block.nx + 1) * (block.ny + 1)];
+		for (int icell = 0; icell < (block.nx + 1) * (block.ny + 1); icell++)
+		{
+			node2d[iblock][icell] = new double[2];
+		}
+	}
+
+
+
+	// section v
+	if (validation_marker == 399)
+	{
+		// read geometries
+		cout << "Geometry section not implemented. Stopping." << endl;
+		exit(-1);
+	}
+	// end of geometrie
+
+	////////////////////////
+	// end of header section
+	////////////////////////
+	cout << "total " << zone_num << " zones are read in" << endl;
+	cout << "the total cell for all zones is " << sum_mesh_number << endl;
+	// Data section
+	assert(validation_marker == 357.0);
+	for (size_t index = 0; index < zone_num; index++)
+	{
+		in.read((char*)&validation_marker, sizeof(validation_marker));
+		assert(validation_marker == 299.0);
+		int format[3];
+		for (int i_var = 0; i_var < number_variables; i_var++)
+		{
+			in.read((char*)&format[i_var], sizeof(format[i_var]));
+			if (format[i_var] != 2)
+			{
+				cout << "we only support double pricision coordinates now" << endl;
+				exit(0);
+			}
+		}
+		int has_passive_variables;
+		in.read((char*)&has_passive_variables, sizeof(has_passive_variables));
+
+		if (has_passive_variables)
+		{
+			cout << "we donnt support has_passive_variables != 0 now" << endl;
+			exit(0);
+		}
+		int has_variable_sharing;
+		in.read((char*)&has_variable_sharing, sizeof(has_variable_sharing));
+
+		if (has_variable_sharing)
+		{
+			cout << "we donnt support has_variable_sharing != 0 now" << endl;
+			exit(0);
+		}
+		int zone_share_connectivity;
+		in.read((char*)&zone_share_connectivity, sizeof(zone_share_connectivity));
+
+		double  min_value[3], max_value[3];
+
+		for (int i = 0; i < number_variables; i++)
+		{
+			in.read((char*)&min_value[i], sizeof(min_value[i]));
+			in.read((char*)&max_value[i], sizeof(max_value[i]));
+		}
+
+		for (int i_var = 0; i_var < number_variables; i_var++)
+		{
+			for (int k = 0; k < kmax[index]; k++)
+			{
+				for (int j = 0; j < jmax[index]; j++)
+				{
+					for (int i = 0; i < imax[index]; i++)
+					{
+						double coord;
+						in.read((char*)&coord, sizeof(coord));
+						int icell = block.ghost + i;
+						int jcell = block.ghost + j;
+						int inode = icell * (block.ny + 1) + jcell;
+						if (i_var < 2)
+						{
+							node2d[index][inode][i_var] = coord;
+						}
+
+					}
+				}
+			}
+		}
+
+	}
+
+	create_mirror_ghost_cell(block, node2d[0]);
+
+	//once we have the dimension infomation, we can use it 
+	// to allocate memory
+	*fluids = Setfluid(block);
+	// then the interfaces reconstructioned value, (N+1)*(N+1)
+	*xinterfaces = Setinterface_array(block);
+	*yinterfaces = Setinterface_array(block);
+	// then the flux, which the number is identical to interfaces
+	*xfluxes = Setflux_gauss_array(block);
+	*yfluxes = Setflux_gauss_array(block);
+
+	//we shall use the node coordinate to compute the geometric infos for cells/interfaces
+
+
+	SetNonUniformMesh(block, fluids[0],
+		xinterfaces[0], yinterfaces[0],
+		xfluxes[0], yfluxes[0], node2d[0]);
+
+}
+
+void read_ascii_to_string(string& value, ifstream& in)
+{
+	int ascii(-1);
+
+	in.read((char*)&ascii, sizeof(ascii));
+	while (ascii != 0)
+	{
+		char temp = (char)ascii;
+		value.append(sizeof(char), temp);
+		in.read(reinterpret_cast<char*>(&ascii), sizeof(ascii));
+	}
+}
+
+void create_mirror_ghost_cell(Block2d& block, double** node2d)
+{
+
+	//create the coordinate for ghost nodes
+	for (int j = 0; j < block.ny + 1; j++)
+	{
+		for (int i = 0; i < block.nx + 1; i++)
+		{
+			int inode = i * (block.ny + 1) + j;
+			int base_node = 0, bc_node = 0;
+			int base_index_x = i, base_index_y = j;
+			int bc_index_x = i, bc_index_y = j;
+			if (i < block.ghost)
+			{
+				base_index_x = 2 * block.ghost - i;
+				bc_index_x = block.ghost;
+			}
+			else if (i > block.nx - block.ghost)
+			{
+				base_index_x = 2 * (block.nx - block.ghost) - i;
+				bc_index_x = block.nx - block.ghost;
+			}
+			if (j < block.ghost)
+			{
+				base_index_y = 2 * block.ghost - j;
+				bc_index_y = block.ghost;
+			}
+			else if (j > block.ny - block.ghost)
+			{
+				base_index_y = 2 * (block.ny - block.ghost) - j;
+				bc_index_y = block.ny - block.ghost;
+			}
+			if (base_index_x != i || base_index_y != j)
+			{
+				base_node = base_index_x * (block.ny + 1) + base_index_y;
+				bc_node = bc_index_x * (block.ny + 1) + bc_index_y;
+
+				//cacluate mirror point
+				node2d[inode][0] = 2 * node2d[bc_node][0] - node2d[base_node][0];
+				node2d[inode][1] = 2 * node2d[bc_node][1] - node2d[base_node][1];
+			}
+		}
+	}
+}

code/mesh_read.h

@@ -0,0 +1,10 @@
+#pragma once
+#include"fvm_gks2d.h"
+#include"output.h"
+using namespace std;
+string add_mesh_directory_modify_for_linux();
+void Read_structured_mesh(string& filename,
+	Fluid2d** fluids, Interface2d** xinterfaces, Interface2d** yinterfaces,
+	Flux2d_gauss*** xfluxes, Flux2d_gauss*** yfluxes, Block2d& block);
+void read_ascii_to_string(string& value, ifstream& in);
+void create_mirror_ghost_cell(Block2d& block, double** node2d);

code/output.cpp

@@ -0,0 +1,823 @@
+#include"output.h"
+
+
+
+void make_directory_for_result()
+{
+	string newfolder = "result";
+
+#if defined(_WIN32)
+	//一个简单的宏定义判断，如果是win系统 则运行以下代码
+	if (_mkdir(newfolder.c_str()) != 0)
+	{
+		_mkdir(newfolder.c_str());
+	}
+	else
+	{
+		cout << "fail to create new folder for result file" << endl;
+	}
+	//end
+#else 
+	//如果是linux系统 则运行以下代码
+	mkdir(newfolder.c_str(), 0777);
+	if (mkdir(newfolder.c_str(), 0777) != 0)
+	{
+		mkdir(newfolder.c_str(), 0777);
+	}
+	else
+	{
+		cout << "fail to create new folder for result file" << endl;
+	}
+	//end
+#endif
+
+}
+
+void couthere()
+{
+	// to test that the program can run up to here
+	cout << "here" << endl;
+}
+void cinhere()
+{
+	// to test taht the program can run up here
+	//and only show once
+	double a;
+	cout << "the program run here" << endl;
+	cin >> a;
+}
+
+void output_error_form(double CFL, double dt_ratio, int mesh_set, int* mesh_number, double** error)
+{
+	cout << "Accuracy-residual-file-output" << endl;
+
+	ofstream error_out("result/error.txt");
+	error_out << "the CFL number is " << CFL << endl;
+	error_out << "the dt ratio over dx is " << dt_ratio << endl;
+
+
+	double** order = new double* [mesh_set - 1];
+	for (int i = 0; i < mesh_set - 1; i++)
+	{
+		order[i] = new double[3];
+		for (int j = 0; j < 3; j++)
+		{
+			order[i][j] = log(error[i][j] / error[i + 1][j]) / log(2);
+		}
+	}
+
+
+	for (int i = 0; i < mesh_set; i++)
+	{
+		if (i == 0)
+		{
+			error_out << "1/" << mesh_number[i]
+				<< scientific << setprecision(6)
+				<< " & " << error[i][0] << " & ~"
+				<< " & " << error[i][1] << " & ~"
+				<< " & " << error[i][2] << " & ~"
+				<< " \\\\ " << endl;
+		}
+		else
+		{
+			error_out << "1/" << mesh_number[i]
+				<< " & " << scientific << setprecision(6) << error[i][0];
+			error_out << " & " << fixed << setprecision(2) << order[i - 1][0];
+			error_out << " & " << scientific << setprecision(6) << error[i][1];
+			error_out << " & " << fixed << setprecision(2) << order[i - 1][1];
+			error_out << " & " << scientific << setprecision(6) << error[i][2];
+			error_out << " & " << fixed << setprecision(2) << order[i - 1][2];
+			error_out << " \\\\ " << endl;
+		}
+
+	}
+	error_out.close();
+}
+
+
+void output1d(Fluid1d* fluids, Block1d block)
+{
+	stringstream name;
+	name << "result/Result1D-" << setfill('0') << setw(5) << block.step << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << "variables = x,density,u,pressure,temperature,entropy,Ma" << endl;
+	out << "zone i = " << block.nodex << ", F=POINT" << endl;
+
+	//before output data, converting conservative varabile to primitive variable
+	Convar_to_primvar(fluids, block);
+	//output the data
+	for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+	{
+		out << fluids[i].cx << " "
+			<< setprecision(10)
+			<< fluids[i].primvar[0] << " "
+			<< fluids[i].primvar[1] << " "
+			<< fluids[i].primvar[2] << " "
+			<< Temperature(fluids[i].primvar[0], fluids[i].primvar[2]) << " "
+			<< entropy(fluids[i].primvar[0], fluids[i].primvar[2]) << " "
+			<< sqrt(pow(fluids[i].primvar[1], 2)) / Soundspeed(fluids[i].primvar[0], fluids[i].primvar[2]) << " "
+			<< endl;
+	}
+	out.close();
+}
+
+void continue_compute(Fluid2d* fluids, Block2d& block)
+{
+	fstream in("result/continue-in.txt", ios::in);
+	if (!in.is_open())
+	{
+		cout << "cannot find continue.txt" << endl;
+		cout << "a new case will start" << endl;
+		return;
+	}
+	int in_mesh[2];
+	in >> in_mesh[0] >> in_mesh[1];
+
+	if (in_mesh[0] == block.nodex && in_mesh[1] == block.nodey)
+	{
+
+		fstream dat("result/continue-in.plt", ios::in);
+		if (!dat.is_open())
+		{
+			cout << "cannot find continue.plt, which contains fluid data" << endl;
+			cout << "a new case will start" << endl;
+			return;
+		}
+		string s;
+		getline(dat, s);
+		dat.ignore(256, '\n');
+
+		if (block.ghost <= 0 || block.ghost >= 5)
+		{
+			cout << "a wrong ghost set up initializing" << endl;
+		}
+
+		//cell avg information
+		for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+		{
+			for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+			{
+				int index = i * block.ny + j;
+				dat >> fluids[index].coordx
+					>> fluids[index].coordy;
+
+				for (int var = 0; var < 4; var++)
+				{
+					dat >> fluids[index].convar[var];
+				}
+
+
+
+			}
+
+		}
+
+
+		in >> block.t;
+		cout << "successfully read in continue-compute-file in t= " << block.t << endl;
+
+		Convar_to_Primvar(fluids, block);
+		cout << "change convar to primvar" << endl;
+	}
+	else
+	{
+		cout << "mesh size is not same" << endl;
+		cout << "a new case will start" << endl;
+		return;
+	}
+	in.close();
+	return;
+}
+void output2d(Fluid2d* fluids, Block2d block)
+{
+	//for create a output file, pls include the follow standard lib
+	//#include <fstream>
+	//#include <sstream>
+	//#include<iomanip>
+	output_record(block);
+	stringstream name;
+	name << "result/Result2D-" << setfill('0') << setw(5) << block.step << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << "# CFL number is " << block.CFL << endl;
+	out << "# tau type (0 refer euler, 1 refer NS) is " << tau_type << endl;
+	if (tau_type == Euler)
+	{
+		out << "# c1_euler " << c1_euler << " c2_euler " << c2_euler << endl;
+	}
+	if (tau_type == NS)
+	{
+		out << "# Mu " << Mu << " Nu " << Nu << endl;
+	}
+	out << "# time marching strategy is " << block.stages << " stage method" << endl;
+	out << "variables = x,y,density,u,v,pressure,temperature,entropy,Ma" << endl;
+	out << "zone i = " << block.nodex << " " << ", j = " << block.nodey << ", F=POINT" << endl;
+
+	//before output data, converting conservative varabile to primitive variable
+	Convar_to_primvar(fluids, block);
+	//output the data
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+
+			out << fluids[i * block.ny + j].coordx << " "
+				<< fluids[i * block.ny + j].coordy << " "
+				<< setprecision(10)
+				<< fluids[i * block.ny + j].primvar[0] << " "
+				<< fluids[i * block.ny + j].primvar[1] << " "
+				<< fluids[i * block.ny + j].primvar[2] << " "
+				<< fluids[i * block.ny + j].primvar[3] << " "
+				<< Temperature(fluids[i * block.ny + j].primvar[0], fluids[i * block.ny + j].primvar[3]) << " "
+				<< entropy(fluids[i * block.ny + j].primvar[0], fluids[i * block.ny + j].primvar[3]) << " "
+				<< sqrt(pow(fluids[i * block.ny + j].primvar[1], 2) + pow(fluids[i * block.ny + j].primvar[2], 2)) /
+				Soundspeed(fluids[i * block.ny + j].primvar[0], fluids[i * block.ny + j].primvar[3]) << " "
+				<< endl;
+		}
+	}
+
+	//fclose(fp);//关闭文件
+	out.close();
+
+}
+void output2d_center(Fluid2d* fluids, Block2d block)
+{
+	output_record(block);
+	stringstream name;
+	name << "result/result2D-" << setfill('0') << setw(5) << block.step << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+
+
+	out << "# CFL number is " << block.CFL << endl;
+	out << "# tau type (0 refer euler, 1 refer NS) is " << tau_type << endl;
+	if (tau_type == Euler)
+	{
+		out << "# c1_euler " << c1_euler << " c2_euler " << c2_euler << endl;
+	}
+	if (tau_type == NS)
+	{
+		out << "# Mu " << Mu << " Nu " << Nu << endl;
+	}
+	out << "# time marching strategy is " << block.stages << " stage method" << endl;
+	out << "variables = x,y,density,u,v,pressure,temperature,entropy,Ma" << endl;
+	out << "zone i = " << block.nodex + 1 << " " << ", j = " << block.nodey + 1 << ", DATAPACKING=BLOCK, VARLOCATION=([3-9]=CELLCENTERED)" << endl;
+
+	//output locations
+	for (int j = block.ghost; j < block.ghost + block.nodey + 1; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex + 1; i++)
+		{
+			out << fluids[i * block.ny + j].node[0] << " ";
+		}
+		out << endl;
+	}
+
+	for (int j = block.ghost; j < block.ghost + block.nodey + 1; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex + 1; i++)
+		{
+			out << fluids[i * block.ny + j].node[1] << " ";
+		}
+		out << endl;
+	}
+
+
+	//before output data, converting conservative varabile to primitive variable
+	Convar_to_primvar(fluids, block);
+	//output the data
+
+	//output four primary variables
+	for (int var = 0; var < 4; var++)
+	{
+		for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+		{
+			for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+			{
+				out << setprecision(10)
+					<< fluids[i * block.ny + j].primvar[var] << " ";
+			}
+			out << endl;
+		}
+	}
+	//output temperature
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+			out << setprecision(10)
+				<< Temperature(fluids[i * block.ny + j].primvar[0], fluids[i * block.ny + j].primvar[3]) << " ";
+		}
+		out << endl;
+	}
+	//output entropy
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+			out << setprecision(10)
+				<< entropy(fluids[i * block.ny + j].primvar[0], fluids[i * block.ny + j].primvar[3]) << " ";
+		}
+		out << endl;
+	}
+
+	//output Ma number
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+			out << setprecision(10)
+				<< sqrt(pow(fluids[i * block.ny + j].primvar[1], 2))
+				/ Soundspeed(fluids[i * block.ny + j].primvar[0], fluids[i * block.ny + j].primvar[3]) << " ";
+		}
+		out << endl;
+	}
+
+	//fclose(fp);//关闭文件
+	out.close();
+}
+
+void output2d_blow_up(Fluid2d* fluids, Block2d block)
+{
+	output_record(block);
+	stringstream name;
+	name << "result/BlowUp-2D-" << setfill('0') << setw(5) << block.step << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << "# CFL number is " << block.CFL << endl;
+	out << "# tau type (0 refer euler, 1 refer NS) is " << tau_type << endl;
+	if (tau_type == Euler)
+	{
+		out << "# c1_euler " << c1_euler << " c2_euler " << c2_euler << endl;
+	}
+	if (tau_type == NS)
+	{
+		out << "# Mu " << Mu << " Nu " << Nu << endl;
+	}
+	out << "# time marching strategy is " << block.stages << " stage method" << endl;
+	out << "variables = x,y,blow_point,density,u,v,pressure" << endl;
+	out << "zone i = " << block.nodex + 1 << " " << ", j = " << block.nodey + 1 << ", DATAPACKING=BLOCK, VARLOCATION=([3-7]=CELLCENTERED)" << endl;
+
+	//output locations
+	for (int j = block.ghost; j < block.ghost + block.nodey + 1; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex + 1; i++)
+		{
+			out << fluids[i * block.ny + j].node[0] << " ";
+		}
+		out << endl;
+	}
+
+	for (int j = block.ghost; j < block.ghost + block.nodey + 1; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex + 1; i++)
+		{
+			out << fluids[i * block.ny + j].node[1] << " ";
+		}
+		out << endl;
+	}
+
+
+	//before output data, converting conservative varabile to primitive variable
+	Convar_to_primvar(fluids, block);
+	//output the data
+
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+			if (negative_density_or_pressure(fluids[i * block.ny + j].primvar))
+			{
+				out << 1.0 << " ";
+			}
+			else
+			{
+				out << 0.0 << " ";
+			}
+		}
+		out << endl;
+	}
+	//output four primary variables
+	for (int var = 0; var < 4; var++)
+	{
+		for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+		{
+			for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+			{
+				if (negative_density_or_pressure(fluids[i * block.ny + j].primvar))
+				{
+					out << -1.0 << " ";
+				}
+				else
+				{
+					out << setprecision(10)
+						<< fluids[i * block.ny + j].primvar[var] << " ";
+				}
+			}
+			out << endl;
+		}
+	}
+
+
+	//fclose(fp);//关闭文件
+	out.close();
+}
+void output_record(Block2d block)
+{
+	//输出一下二维的模拟工况，为了方便记录之用
+	stringstream name;
+	name << "result/Record2D-" << setfill('0') << setw(5) << block.step << ".txt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << "# This is a two dimensional result" << endl;
+	if (block.uniform == true)
+	{
+		out << "# uniform mesh is used for calculation" << endl;
+		out << "# left region coordinate " << endl;
+		out << block.left << endl;
+		out << "# right region coordinate " << endl;
+		out << block.right << endl;
+		out << "# down region coordinate " << endl;
+		out << block.down << endl;
+		out << "# up region coordinate " << endl;
+		out << block.up << endl;
+	}
+	else
+	{
+		out << "# non-uniform mesh is used for calculation" << endl;
+	}
+	out << "# ghost cell number is " << endl;
+	out << block.ghost << endl;
+
+	out << "# Mesh size nodex and nodey" << endl;
+	out << block.nodex << " " << block.nodey << endl;
+	out << "# CFL number is " << endl;
+	out << block.CFL << endl;
+	out << "# simulation time is " << endl;
+	out << block.t << endl;
+	out << "# simulation step is " << endl;
+	out << block.step << endl;
+	out << "# tau type (0 refer euler, 1 refer NS) is " << endl;
+	out << tau_type << endl;
+
+	if (tau_type == Euler)
+	{
+		out << "# c1_euler " << " c2_euler " << endl;
+		out << c1_euler << " " << c2_euler << endl;
+	}
+
+	if (tau_type == NS)
+	{
+		out << "# Mu " << " Nu " << endl;
+		out << Mu << " " << Nu << endl;
+	}
+	out << "# the number of gaussian point it use: " << endl;
+	out << gausspoint << endl;
+
+	out << "# normal_reconstruction_for_left_right_states is ";
+	if (cellreconstruction_2D_normal == First_order_normal)
+	{
+		out << "FirstOrder_splitting " << endl;
+	}
+	else if (cellreconstruction_2D_normal == Vanleer_normal)
+	{
+		out << "Vanleer_splitting " << endl;
+	}
+	else
+	{
+		out << "not_specified" << endl;
+	}
+
+
+	out << "# normal_reconstruction_for_g0_states is ";
+	if (g0reconstruction_2D_normal == Center_collision_normal)
+	{
+		out << "Center_collision_normal " << endl;
+	}
+	else if (g0reconstruction_2D_normal == Center_3rd_normal)
+	{
+		out << "Center_3rd_normal " << endl;
+	}
+	else
+	{
+		out << "not_specified" << endl;
+	}
+
+	out << "# tangential_reconstruction_for_left_right_states is ";
+	if (cellreconstruction_2D_tangent == Vanleer_tangent)
+	{
+		out << "Vanleer_tangent " << endl;
+	}
+	else if (cellreconstruction_2D_tangent == First_order_tangent)
+	{
+		out << "First_order_tangent " << endl;
+	}
+	else
+	{
+		out << "not_specified" << endl;
+	}
+
+	out << "# g0reconstruction_2D_multi is ";
+	if (g0reconstruction_2D_tangent == Center_all_collision_multi)
+	{
+		out << "Center_all_collision_tangent " << endl;
+	}
+	else if (g0reconstruction_2D_tangent == Center_3rd_tangent)
+	{
+		out << "Center_all_collision_tangent " << endl;
+	}
+	else
+	{
+		out << "not_specified" << endl;
+	}
+
+	out << "# time marching strategy is ";
+	if (timecoe_list_2d == S1O1_2D)
+	{
+		out << "S1O1 " << endl;
+	}
+	if (timecoe_list_2d == S2O4_2D)
+	{
+		out << "S2O4 " << endl;
+	}
+	if (timecoe_list_2d == RK4_2D)
+	{
+		out << "RK4 " << endl;
+	}
+
+	out << "# the flux function is ";
+	if (flux_function_2d == GKS2D_smooth)
+	{
+		out << "gks2nd_smooth" << endl;
+	}
+	if (flux_function_2d == GKS2D)
+	{
+		if (gks2dsolver == kfvs1st_2d)
+		{
+			out << "kfvs1st " << endl;
+		}
+		if (gks2dsolver == kfvs2nd_2d)
+		{
+			out << "kfvs2nd " << endl;
+		}
+		if (gks2dsolver == gks1st_2d)
+		{
+			out << "gks1st " << endl;
+		}
+		if (gks2dsolver == gks2nd_2d)
+		{
+			out << "gks2nd " << endl;
+		}
+	}
+
+	out << "# the limitation for current scheme is " << endl;
+	out.close();
+}
+void output2d_binary(Fluid2d* fluids, Block2d block)
+{
+	cout << "Output binary result..." << endl;
+	//before output data, converting conservative varabile to primitive variable
+	Convar_to_primvar(fluids, block);
+
+	output_record(block);
+	stringstream name;
+	name << "result/Result2D-" << setfill('0') << setw(5) << block.step << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out;
+	out.open(s, ios::out | ios::binary);
+
+	int IMax = block.nodex;
+	int JMax = block.nodey;
+	int KMax = 1;
+	int NumVar = 7;
+
+	char Title[] = "structure2d";
+
+
+	char Varname1[] = "X";
+	char Varname2[] = "Y";
+	char Varname3[] = "Z";
+	char Varname4[] = "density";
+	char Varname5[] = "u";
+	char Varname6[] = "v";
+	char Varname7[] = "pressure";
+
+	char Zonename1[] = "Zone 001";
+
+	float ZONEMARKER = 299.0;
+
+	float EOHMARKER = 357.0;
+
+	//==============Header Secontion =================//
+
+	//------1.1 Magic number, Version number
+
+	char MagicNumber[] = "#!TDV101";
+
+	// the magic number must be output with 8 lengths without NULL end
+	out.write(MagicNumber, 8);
+	//---- - 1.2.Integer value of 1.----------------------------------------------------------
+
+	int IntegerValue = 1;
+	out.write((char*)&IntegerValue, sizeof(IntegerValue));
+
+	//---- - 1.3.Title and variable names.------------------------------------------------ -
+
+	//---- - 1.3.1.The TITLE.
+	write_character_into_plt(Title, out);
+	//---- - 1.3.2 Number of variables(NumVar) in the c_strfile.
+
+
+
+	out.write((char*)&NumVar, sizeof(NumVar));
+
+
+	//------1.3.3 Variable names.N = L[1] + L[2] + ....L[NumVar]
+
+	write_character_into_plt(Varname1, out);
+	write_character_into_plt(Varname2, out);
+	write_character_into_plt(Varname3, out);
+	write_character_into_plt(Varname4, out);
+	write_character_into_plt(Varname5, out);
+	write_character_into_plt(Varname6, out);
+	write_character_into_plt(Varname7, out);
+
+	//---- - 1.4.Zones------------------------------------------------------------------ -
+
+	//--------Zone marker.Value = 299.0
+
+	out.write((char*)&ZONEMARKER, sizeof(ZONEMARKER));
+
+	//--------Zone name.
+	write_character_into_plt(Zonename1, out);
+	//--------Zone color
+
+	int ZoneColor = -1;
+
+	out.write((char*)&ZoneColor, sizeof(ZoneColor));
+
+	//--------ZoneType
+
+	int ZoneType = 0;
+
+	out.write((char*)&ZoneType, sizeof(ZoneType));
+
+	//--------DaraPacking 0=Block, 1=Point
+
+	int DaraPacking = 1;
+
+	out.write((char*)&DaraPacking, sizeof(DaraPacking));
+	//--------Specify Var Location. 0 = Don't specify, all c_str is located at the nodes. 1 = Specify
+
+	int SpecifyVarLocation = 0;
+
+	out.write((char*)&SpecifyVarLocation, sizeof(SpecifyVarLocation));
+	//--------Number of user defined face neighbor connections(value >= 0)
+
+	int NumOfNeighbor = 0;
+
+	out.write((char*)&NumOfNeighbor, sizeof(NumOfNeighbor));
+
+	//-------- - IMax, JMax, KMax
+
+	out.write((char*)&IMax, sizeof(IMax));
+	out.write((char*)&JMax, sizeof(JMax));
+	out.write((char*)&KMax, sizeof(KMax));
+
+	//----------// -1 = Auxiliary name / value pair to follow 0 = No more Auxiliar name / value pairs.
+
+	int AuxiliaryName = 0;
+
+	out.write((char*)&AuxiliaryName, sizeof(AuxiliaryName));
+
+	//----I HEADER OVER--------------------------------------------------------------------------------------------
+
+	//=============================Geometries section=======================
+
+	//=============================Text section======================
+
+	// EOHMARKER, value = 357.0
+
+	out.write((char*)&EOHMARKER, sizeof(EOHMARKER));
+
+	//================II.Data section===============//
+
+	//---- - 2.1 zone---------------------------------------------------------------------- -
+
+	out.write((char*)&ZONEMARKER, sizeof(ZONEMARKER));
+
+	//--------variable c_str format, 1 = Float, 2 = Double, 3 = LongInt, 4 = ShortInt, 5 = Byte, 6 = Bit
+
+	int fomat[20];
+	for (int ifomat = 0; ifomat < 20; ifomat++)
+	{
+		fomat[ifomat] = 2;
+	}
+	for (int ifomat = 0; ifomat < NumVar; ifomat++)
+	{
+		out.write((char*)&fomat[ifomat], sizeof(fomat[ifomat]));
+	}
+
+
+	//--------Has variable sharing 0 = no, 1 = yes.
+
+	int HasVarSharing = 0;
+
+	out.write((char*)&HasVarSharing, sizeof(HasVarSharing));
+
+	//----------Zone number to share connectivity list with(-1 = no sharing).
+
+	int ZoneNumToShareConnectivity = -1;
+
+	out.write((char*)&ZoneNumToShareConnectivity, sizeof(ZoneNumToShareConnectivity));
+	//----------Zone c_str.Each variable is in c_str format asspecified above.
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+			double coordz = 0.0;
+			out.write((char*)&fluids[i * block.ny + j].coordx, sizeof(double));
+			out.write((char*)&fluids[i * block.ny + j].coordy, sizeof(double));
+			out.write((char*)&coordz, sizeof(double));
+			out.write((char*)&fluids[i * block.ny + j].primvar[0], sizeof(double));
+			out.write((char*)&fluids[i * block.ny + j].primvar[1], sizeof(double));
+			out.write((char*)&fluids[i * block.ny + j].primvar[2], sizeof(double));
+			out.write((char*)&fluids[i * block.ny + j].primvar[3], sizeof(double));
+		}
+
+	}
+	out.close();
+
+}
+void write_character_into_plt(char* str, ofstream& out)
+{
+	int value = 0;
+	while ((*str) != '\0')
+	{
+		value = (int)*str;
+		out.write((char*)&value, sizeof(value));
+		str++;
+	}
+
+	char null_char[] = "";
+	value = (int)*null_char;
+	out.write((char*)&value, sizeof(value));
+}
+void continue_compute_output(Fluid2d* fluids, Block2d block)
+{
+	//for create a output file, pls include the follow standard lib
+	//#include <fstream>
+	//#include <sstream>
+	//#include<iomanip>
+	continue_compute_output_param(block);
+	stringstream name;
+	name << "result/continue-out" << ".plt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << "variables = x,y,convar0,convar1,convar2,convar3";
+	out << endl;
+
+	out << "zone i = " << block.nodex << " " << ", j = " << block.nodey << ", F=POINT" << endl;
+
+	//before output data, converting conservative varabile to primitive variable
+	Convar_to_primvar(fluids, block);
+
+	//output the data
+	for (int j = block.ghost; j < block.ghost + block.nodey; j++)
+	{
+		for (int i = block.ghost; i < block.ghost + block.nodex; i++)
+		{
+			out << fluids[i * block.ny + j].coordx << " "
+				<< fluids[i * block.ny + j].coordy << " "
+				<< setprecision(10);
+			for (int var = 0; var < 4; var++)
+			{
+				out << fluids[i * block.ny + j].convar[var] << " ";
+			}
+
+			out << endl;
+		}
+	}
+
+	out.close();
+}
+void continue_compute_output_param(Block2d block)
+{
+	stringstream name;
+	name << "result/continue-out" << ".txt" << endl;
+	string s;
+	name >> s;
+	ofstream out(s);
+	out << block.nodex << " " << block.nodey << endl; //mesh point
+	out << setprecision(14);
+	out << block.t << endl;//simulation time
+}
+

code/output.h

@@ -0,0 +1,49 @@
+#pragma once
+#include"fvm_gks2d.h"
+#include"time.h"
+#include<string.h>
+#include<iostream>
+#include<fstream>
+#include<sstream>
+#include<iomanip>
+#if defined(_WIN32)
+#include<direct.h>
+#else 
+#include<sys/stat.h>
+//#include<sys/types.h>
+#endif
+struct Runtime {
+	//for record wall clock time
+	clock_t start_initial;
+	clock_t finish_initial;
+	clock_t start_compute;
+	clock_t finish_compute;
+	Runtime()
+	{
+		memset(this, 0, sizeof(Runtime));
+	}
+};
+
+void make_directory_for_result();
+
+void couthere();
+void cinhere();
+
+void output_error_form(double CFL, double dt_ratio, int mesh_set, int* mesh_number, double** error);
+
+//tecplot 一维输出
+void output1d(Fluid1d* fluids, Block1d block);
+void continue_compute(Fluid2d* fluids, Block2d& block);
+//tecplot 一维输出，单元中心点值
+void output2d(Fluid2d* fluids, Block2d block);
+//tecplot 一维输出，单元值
+void output2d_center(Fluid2d* fluids, Block2d block);
+//输出哪里爆掉了
+void output2d_blow_up(Fluid2d* fluids, Block2d block);
+void output_record(Block2d block);
+void output2d_binary(Fluid2d* fluids, Block2d block);
+void write_character_into_plt(char* str, ofstream& out);
+void continue_compute_output(Fluid2d* fluids, Block2d block);
+void continue_compute_output_param(Block2d block);
+
+

code/riemann_problem.cpp

@@ -0,0 +1,448 @@
+#include"riemann_problem.h"
+
+void riemann_problem_1d()
+{
+	Runtime runtime;//statement for recording the running time
+	runtime.start_initial = clock();
+
+	Block1d block;
+	block.nodex = 100;
+	block.ghost = 4;
+
+	double tstop = 0.2;
+	block.CFL = 0.5;
+	Fluid1d* bcvalue = new Fluid1d[2];
+	K = 4;
+	Gamma = 1.4;
+	R_gas = 1.0;
+
+	gks1dsolver = gks2nd;
+	//g0type = collisionn;
+	tau_type = Euler;
+	//Smooth = false;
+	c1_euler = 0.05;
+	c2_euler = 1;
+
+	//prepare the boundary condtion function
+	BoundaryCondition leftboundary(0);
+	BoundaryCondition rightboundary(0);
+	leftboundary = free_boundary_left;
+	rightboundary = free_boundary_right;
+	//prepare the reconstruction function
+
+	cellreconstruction = Vanleer;
+	wenotype = wenoz;
+	reconstruction_variable = characteristic;
+	g0reconstruction = Center_3rd;
+	is_reduce_order_warning = true;
+	//prepare the flux function
+	flux_function = GKS;
+	//prepare time marching stratedgy
+	timecoe_list = S1O2;
+	Initial_stages(block);
+
+	// allocate memory for 1-D fluid field
+	// in a standard finite element method, we have 
+	// first the cell average value, N
+	block.nx = block.nodex + 2 * block.ghost;
+	block.nodex_begin = block.ghost;
+	block.nodex_end = block.nodex + block.ghost - 1;
+	Fluid1d* fluids = new Fluid1d[block.nx];
+	// then the interfaces reconstructioned value, N+1
+	Interface1d* interfaces = new Interface1d[block.nx + 1];
+	// then the flux, which the number is identical to interfaces
+	Flux1d** fluxes = Setflux_array(block);
+	//end the allocate memory part
+
+	//bulid or read mesh,
+	//since the mesh is all structured from left and right,
+	//there is no need for buliding the complex topology between cells and interfaces
+	//just using the index for address searching
+
+	block.left = 0.0;
+	block.right = 1.0;
+	block.dx = (block.right - block.left) / block.nodex;
+	//set the uniform geometry information
+	SetUniformMesh(block, fluids, interfaces, fluxes);
+
+	//ended mesh part
+
+	// then its about initializing, lets first initialize a sod test case
+	//you can initialize whatever kind of 1d test case as you like
+	Fluid1d zone1; zone1.primvar[0] = 1.0; zone1.primvar[1] = 0.0; zone1.primvar[2] = 1.0;
+	Fluid1d zone2; zone2.primvar[0] = 0.125; zone2.primvar[1] = 0.0; zone2.primvar[2] = 0.1;
+
+	//Fluid1d zone1; zone1.primvar[0] = 1.0; zone1.primvar[1] = -2.0; zone1.primvar[2] = 0.4;
+	//Fluid1d zone2; zone2.primvar[0] = 1.0; zone2.primvar[1] = 2.0; zone2.primvar[2] = 0.4;
+	ICfor1dRM(fluids, zone1, zone2, block);
+	//initializing part end
+
+
+	//then we are about to do the loop
+	block.t = 0;//the current simulation time
+	block.step = 0; //the current step
+
+	int inputstep = 1;//input a certain step,
+					  //initialize inputstep=1, to avoid a 0 result
+	runtime.finish_initial = clock();
+	while (block.t < tstop)
+	{
+		// assume you are using command window,
+		// you can specify a running step conveniently
+		if (block.step % inputstep == 0)
+		{
+			cout << "pls cin interation step, if input is 0, then the program will exit " << endl;
+			cin >> inputstep;
+			if (inputstep == 0)
+			{
+				output1d(fluids, block);
+				break;
+			}
+		}
+		if (runtime.start_compute == 0.0)
+		{
+			runtime.start_compute = clock();
+			cout << "runtime-start " << endl;
+		}
+
+		//Copy the fluid vales to fluid old
+		CopyFluid_new_to_old(fluids, block);
+		//determine the cfl condtion
+		block.dt = Get_CFL(block, fluids, tstop);
+
+
+		for (int i = 0; i < block.stages; i++)
+		{
+			//after determine the cfl condition, let's implement boundary condtion
+			leftboundary(fluids, block, bcvalue[0]);
+			rightboundary(fluids, block, bcvalue[1]);
+			// here the boudary type, you shall go above the search the key words"BoundaryCondition leftboundary;"
+			// to see the pointer to the corresponding function
+
+			//then is reconstruction part, which we separate the left or right reconstrction
+			//and the center reconstruction
+			Reconstruction_within_cell(interfaces, fluids, block);
+			Reconstruction_forg0(interfaces, fluids, block);
+
+			//then is solver part
+			Calculate_flux(fluxes, interfaces, block, i);
+
+			//then is update flux part
+			Update(fluids, fluxes, block, i);
+			//output1d_checking(fluids, interfaces, fluxes, block);
+		}
+
+		block.step++;
+		block.t = block.t + block.dt;
+		if (block.step % 10 == 0)
+		{
+			cout << "step 10 time is " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+		}
+		if ((block.t - tstop) > 0)
+		{
+			runtime.finish_compute = clock();
+			cout << "initializiation time is " << (float)(runtime.finish_initial - runtime.start_initial) / CLOCKS_PER_SEC << "seconds" << endl;
+			cout << "computational time is " << (float)(runtime.finish_compute - runtime.start_compute) / CLOCKS_PER_SEC << "seconds" << endl;
+			output1d(fluids, block);
+			//output1d_checking(fluids, interfaces, fluxes, block);
+		}
+	}
+
+}
+
+void ICfor1dRM(Fluid1d* fluids, Fluid1d zone1, Fluid1d zone2, Block1d block)
+{
+	for (int i = 0; i < block.nx; i++)
+	{
+		if (i < 0.5 * block.nx)
+		{
+			Copy_Array(fluids[i].primvar, zone1.primvar, 3);
+		}
+		else
+		{
+			Copy_Array(fluids[i].primvar, zone2.primvar, 3);
+		}
+
+	}
+
+	for (int i = 0; i < block.nx; i++)
+	{
+		Primvar_to_convar_1D(fluids[i].convar, fluids[i].primvar);
+	}
+}
+
+void riemann_problem_2d()
+{
+
+	Runtime runtime;
+	runtime.start_initial = clock();
+	Block2d block;
+	block.uniform = true;
+	block.nodex = 200;
+	block.nodey = 200;
+	block.ghost = 3;
+
+
+
+	block.CFL = 0.5;
+	Fluid2d* bcvalue = new Fluid2d[4];
+
+	K = 3;
+	Gamma = 1.4;
+
+
+
+
+	//prepare the boundary condtion function
+	BoundaryCondition2dnew leftboundary(0);
+	BoundaryCondition2dnew rightboundary(0);
+	BoundaryCondition2dnew downboundary(0);
+	BoundaryCondition2dnew upboundary(0);
+
+	leftboundary = free_boundary_left;
+	rightboundary = free_boundary_right;
+	downboundary = free_boundary_down;
+	upboundary = free_boundary_up;
+
+	//prepare the reconstruction function
+
+	gausspoint = 1;
+	SetGuassPoint();
+
+	reconstruction_variable = characteristic;
+	wenotype = wenoz;
+
+	cellreconstruction_2D_normal = Vanleer_normal;
+	cellreconstruction_2D_tangent = Vanleer_tangent;
+	g0reconstruction_2D_normal = Center_do_nothing_normal;
+	g0reconstruction_2D_tangent = Center_all_collision_multi;
+
+	is_reduce_order_warning = true; //重构出负后 输出出负的单元
+
+	//prepare the flux function
+	gks2dsolver = gks2nd_2d;
+	tau_type = Euler;
+	c1_euler = 0.05;
+	c2_euler = 1;
+	flux_function_2d = GKS2D;
+
+	//prepare time marching stratedgy
+
+
+	//time coe list must be 2d
+	timecoe_list_2d = S1O2_2D;
+	Initial_stages(block);
+
+
+	// allocate memory for 2-D fluid field
+	// in a standard finite element method, we have 
+	// first the cell average value, N*N
+	block.nx = block.nodex + 2 * block.ghost;
+	block.ny = block.nodey + 2 * block.ghost;
+
+	Fluid2d* fluids = Setfluid(block);
+	// then the interfaces reconstructioned value, (N+1)*(N+1)
+	Interface2d* xinterfaces = Setinterface_array(block);
+	Interface2d* yinterfaces = Setinterface_array(block);
+	// then the flux, which the number is identical to interfaces
+	Flux2d_gauss** xfluxes = Setflux_gauss_array(block);
+	Flux2d_gauss** yfluxes = Setflux_gauss_array(block);
+	//end the allocate memory part
+
+	block.left = 0.0;
+	block.right = 1.0;
+	block.down = 0.0;
+	block.up = 1.0;
+	block.dx = (block.right - block.left) / block.nodex;
+	block.dy = (block.up - block.down) / block.nodey;
+	block.overdx = 1 / block.dx;
+	block.overdy = 1 / block.dy;
+	//set the uniform geometry information
+	SetUniformMesh(block, fluids, xinterfaces, yinterfaces, xfluxes, yfluxes);
+	//ended mesh part
+
+	//RM 1 T=0.2
+	double tstop[] = { 0.2 };
+	double zone1[] = { 0.8, 0, 0, 1.0 };
+	double zone2[]{ 1.0, 0, 0.7276, 1.0 };
+	double zone3[]{ 0.5313, 0, 0, 0.4 };
+	double zone4[]{ 1.0, 0.7276, 0, 1.0 };
+
+	////RM config 1 four rarefaction wave T=0.2
+	//double tstop[] = { 0.2 };
+	//Initial zone1{ 0.1072,-0.7259,-1.4045,0.0439 };
+	//Initial zone2{ 0.2579,0,-1.4045,0.15};
+	//Initial zone3{ 1,0,0,1};
+	//Initial zone4{ 0.5197,-0.7259,0,0.4};
+
+	//double tstop[] = { 0.2 };
+	//Initial zone1{ 1.0, 0, 0, 1.0 };
+	//Initial zone4{ 1.0, 0, 0, 1.0 };
+	//Initial zone3{ 0.125, 0, 0, 0.1 };
+	//Initial zone2{ 0.125, 0, 0, 0.1 };
+
+	//double tstop[] = { 0.2 };
+	//Initial zone1{ 1.0, 0, 0, 1.0 };
+	//Initial zone2{ 1.0, 0, 0, 1.0 };
+	//Initial zone3{ 0.125, 0, 0, 0.1 };
+	//Initial zone4{ 0.125, 0, 0, 0.1 };
+
+	////RM 3 T=0.8 with x,y=0.5
+	//double tstop[] = { 0.4, 0.45, 0.5, 0.55,0.6 };
+	//double p0 = 1.0; // this p0 is for adjust the speed of shear layer
+	//Initial zone1{ 1, -0.75, 0.5, p0 };
+	//Initial zone2{ 3.0, -0.75, -0.5, p0 };
+	//Initial zone3{ 1.0, 0.75, -0.5, p0 };
+	//Initial zone4{ 2.0, 0.75, 0.5, p0 };
+
+	////RM 2 T=0.6 with x,y = 0.7
+	//double tstop[] = { 0.3, 0.4, 0.5 };
+	//Initial zone1{ 0.138, 1.206, 1.206, 0.029 };
+	//Initial zone2{ 0.5323, 0, 1.206, 0.3 };
+	//Initial zone3{ 1.5, 0, 0, 1.5 };
+	//Initial zone4{ 0.5323, 1.206, 0, 0.3 };
+
+	////RM 4 T=0.25 with x,y=0.5
+	//double tstop[] = { 0.25,0.3};
+	//double p0 = 1.0; // this p0 is for adjust the speed of shear layer
+	//Initial zone3{ 1, -0.75, -0.5, p0 };
+	//Initial zone4{ 2.0, -0.75, 0.5, p0 };
+	//Initial zone1{ 1.0, 0.75, 0.5, p0 };
+	//Initial zone2{ 3.0, 0.75, -0.5, p0 };
+
+	IC_for_riemann_2d(fluids, zone1, zone2, zone3, zone4, block);
+
+	runtime.finish_initial = clock();
+	block.t = 0;//the current simulation time
+	block.step = 0; //the current step
+	int tstop_dim = sizeof(tstop) / sizeof(double);
+
+	int inputstep = 1;//input a certain step,
+					  //initialize inputstep=1, to avoid a 0 result
+
+	//在若干个模拟时间输出
+	for (int instant = 0; instant < tstop_dim; ++instant)
+	{
+
+		if (inputstep == 0)
+		{
+			break;
+		}
+		while (block.t < tstop[instant])
+		{
+
+			if (block.step % inputstep == 0)
+			{
+				cout << "pls cin interation step, if input is 0, then the program will exit " << endl;
+				cin >> inputstep;
+				if (inputstep == 0)
+				{
+					output2d_binary(fluids, block);
+					break;
+				}
+			}
+			if (runtime.start_compute == 0.0)
+			{
+				runtime.start_compute = clock();
+				cout << "runtime-start " << endl;
+			}
+
+			CopyFluid_new_to_old(fluids, block);
+
+			block.dt = Get_CFL(block, fluids, tstop[instant]);
+
+			for (int i = 0; i < block.stages; i++)
+			{
+				leftboundary(fluids, block, bcvalue[0]);
+				rightboundary(fluids, block, bcvalue[1]);
+				downboundary(fluids, block, bcvalue[2]);
+				upboundary(fluids, block, bcvalue[3]);
+
+				Convar_to_Primvar(fluids, block);
+
+				Reconstruction_within_cell(xinterfaces, yinterfaces, fluids, block);
+
+				Reconstruction_forg0(xinterfaces, yinterfaces, fluids, block);
+
+				Calculate_flux(xfluxes, yfluxes, xinterfaces, yinterfaces, block, i);
+
+				Update(fluids, xfluxes, yfluxes, block, i);
+
+			}
+
+			block.step++;
+			block.t = block.t + block.dt;
+
+			if (block.step % 10 == 0)
+			{
+				cout << "step 10 time is " << (double)(clock() - runtime.start_compute) / CLOCKS_PER_SEC << endl;
+			}
+
+			if ((block.t - tstop[instant]) > 0)
+			{
+				output2d(fluids, block);
+			}
+			if (block.step % 100 == 0)
+			{
+				output2d(fluids, block);
+			}
+			if (fluids[block.nx / 2 * (block.ny) + block.ny / 2].convar[0] != fluids[block.nx / 2 * (block.ny) + block.ny / 2].convar[0])
+			{
+				cout << "the program blows up!" << endl;
+				output2d(fluids, block);
+				break;
+			}
+		}
+	}
+	runtime.finish_compute = clock();
+	cout << "the total run time is " << (double)(runtime.finish_compute - runtime.start_initial) / CLOCKS_PER_SEC << " second !" << endl;
+	cout << "initializing time is" << (double)(runtime.finish_initial - runtime.start_initial) / CLOCKS_PER_SEC << " second !" << endl;
+	cout << "the pure computational time is" << (double)(runtime.finish_compute - runtime.start_compute) / CLOCKS_PER_SEC << " second !" << endl;
+
+}
+
+void IC_for_riemann_2d(Fluid2d* fluid, double* zone1, double* zone2, double* zone3, double* zone4, Block2d block)
+{
+	double discon[2];
+	discon[0] = 0.5;
+	discon[1] = 0.5;
+	for (int i = block.ghost; i < block.nx - block.ghost; i++)
+	{
+		for (int j = block.ghost; j < block.ny - block.ghost; j++)
+		{
+			if (i < discon[0] * block.nx && j < discon[1] * block.ny)
+			{
+				Copy_Array(fluid[i * block.ny + j].primvar, zone1, 4);
+			}
+			else
+			{
+				if (i >= discon[0] * block.nx && j < discon[1] * block.ny)
+				{
+					Copy_Array(fluid[i * block.ny + j].primvar, zone2, 4);
+				}
+				else
+				{
+					if (i >= discon[0] * block.nx && j >= discon[1] * block.ny)
+					{
+						Copy_Array(fluid[i * block.ny + j].primvar, zone3, 4);
+					}
+					else
+					{
+						Copy_Array(fluid[i * block.ny + j].primvar, zone4, 4);
+					}
+				}
+			}
+
+		}
+
+	}
+
+	for (int i = 0; i < block.nx; i++)
+	{
+		for (int j = 0; j < block.ny; j++)
+		{
+			Primvar_to_convar_2D(fluid[i * block.ny + j].convar, fluid[i * block.ny + j].primvar);
+
+		}
+	}
+}

code/riemann_problem.h

@@ -0,0 +1,8 @@
+#pragma once
+#include"fvm_gks2d.h"
+#include"output.h"
+void riemann_problem_1d();
+void ICfor1dRM(Fluid1d* fluids, Fluid1d zone1, Fluid1d zone2, Block1d block);
+void riemann_problem_2d();
+void IC_for_riemann_2d(Fluid2d* fluid, double* zone1, double* zone2, double* zone3, double* zone4, Block2d block);
+