mech.cpp

this a Navier-Stokes equations solver. It support grids contains of multiple connected rectangles. S

{
	prof.start();

	real ret = solve9();
	
	total_mass = 0;
	for(size_t i = 0; i < count; i++){
		grids[i].update_minmax();
		total_mass += grids[i].mass;
		swap(grids[i].array,grids[i].prev_array);
	}

	last_calc_time = prof.get();
	//printf("mech::solve() %.3f sec, mass %f\n",last_calc_time,total_mass);
	return ret;
}

real mech::solve5()
{
	node diff;

	for(size_t i = 0; i < count; i++){
		const node *p = grids[i].prev_array;
		int width = grids[i].width, height = grids[i].height;
		int lastline = width*(height - 1);

		// top line
		for(int x = 0; x < width; x++){
			const node *left, *right;

			if(x > 0)left = p + x - 1;
			else left = get_bound(&grids[i], x - 1, 0);

			if(x < width - 1)right = p + x + 1;
			else right = get_bound(&grids[i], x + 1, 0);

			solve_equation(grids[i].array + x,
				p + x,
				left,
				get_bound(&grids[i],x,-1),
				right,
				p + x + width,
				grids[i].dx,grids[i].dy,dt);
		}

		// bottom line
		for(int x = 0; x < width; x++){
			const node *left, *right;
			if(x > 0)left = p + x - 1 + lastline;
			else left = get_bound(&grids[i], x - 1, height - 1);

			if(x < width - 1)right = p + x + 1 + lastline;
			else right = get_bound(&grids[i], x + 1, height - 1);

			solve_equation(grids[i].array + x + lastline,
				p + x + lastline,
				left,
				p + x - width + lastline,
				right,
				get_bound(&grids[i], x, height),
				grids[i].dx,grids[i].dy,dt);
		}

		// left line
		for(int y = 1; y < height - 1; y++){
			int idx = y*width;

			solve_equation(grids[i].array + idx,
				p + idx,
				get_bound(&grids[i], -1, y),
				p + idx - width,
				p + idx + 1,
				p + idx + width,
				grids[i].dx,grids[i].dy,dt);
		}

		// right line
		for(int y = 1; y < height - 1; y++){
			int idx = (y + 1)*width - 1;
			solve_equation(grids[i].array + idx,
				p + idx,
				p + idx - 1,
				p + idx - width,
				get_bound(&grids[i], width, y),
				p + idx + width,
				grids[i].dx,grids[i].dy,dt);
		}

		// all other points of grid
		for(int y = 1; y < height - 1; y++)for(int x = 1; x < width - 1; x++){
			int idx = x + y*width;
			const node
				*left   = p + idx - 1,
				*top    = p + idx + 1,
				*right  = p + idx - width,
				*bottom = p + idx + width,
				*center = p + idx;

			solve_equation(grids[i].array + idx,center,left,top,right,bottom,grids[i].dx,grids[i].dy,dt);
		}

	}
	return 0;
}

real mech::solve9()
{
	for(size_t i = 0; i < count; i++){
		const node *p = grids[i].prev_array;
		int width = grids[i].width, height = grids[i].height;
		const node *a[9];
		edge9_t edge(a);

		// top border
		for(int x = 0; x < width; x++){
			for(int j = -1; j < 2; j++)for(int k = -1; k < 2; k++){
				if(x + k < 0 || x + k >= width || j < 0 || j >= height)
					edge(k,j) = get_bound(&grids[i],x + k, j);
				else edge(k,j) = p + x + k + j*width;
			}
			solve_equation9(grids[i].array + x,edge,grids[i].dx,grids[i].dy,dt);
		}

		// bottom border
		for(int x = 0; x < width; x++){
			for(int j = -1; j < 2; j++)for(int k = -1; k < 2; k++){
				if(x + k < 0 || x + k >= width || j + height - 1 < 0 || j + height - 1 >= height)
					edge(k,j) = get_bound(&grids[i],x + k, height -1 + j);
				else edge(k,j) = p + x + k + (height - 1 + j)*width;
			}
			solve_equation9(grids[i].array + x + (height - 1)*width,edge,grids[i].dx,grids[i].dy,dt);
		}

		// left border
		for(int y = 1; y < height - 1; y++){
			for(int j = -1; j < 2; j++)for(int k = -1; k < 2; k++){
				if(k < 0 || k >= width || j + y < 0 || j + y >= height)
					edge(k,j) = get_bound(&grids[i],k, y + j);
				else edge(k,j) = p + k + (y + j)*width;
			}
			solve_equation9(grids[i].array + y*width,edge,grids[i].dx,grids[i].dy,dt);
		}

		// right border
		for(int y = 1; y < height - 1; y++){
			for(int j = -1; j < 2; j++)for(int k = -1; k < 2; k++){
				if(width - 1 + k < 0 || width - 1 + k >= width || j + y < 0 || j + y >= height)
					edge(k,j) = get_bound(&grids[i],width - 1 + k, y + j);
				else edge(k,j) = p + width - 1 + k + (y + j)*width;
			}
			solve_equation9(grids[i].array + width - 1 + y*width,edge,grids[i].dx,grids[i].dy,dt);
		}

		for(int y = 1; y < height - 1; y++)for(int x = 1; x < width - 1; x++){
			int idx = x + y*width;

			for(int j = -1; j < 2; j++)for(int k = -1; k < 2; k++)
				edge(k,j) = p + idx + k + j*width;

			solve_equation9(grids[i].array + idx,edge,grids[i].dx,grids[i].dy,dt);
		}
	}
	return 0;
}

int mech::solve_equation(node *result,
						 const node *center,
						 const node *left,
						 const node *top,
						 const node *right,
						 const node *bottom,
						 real dx,real dy,real dt)
{
	// simple thermal diffuse equation
	result->T = center->T + dt*((right->T - 2*center->T + left->T)/dx/dx + (bottom->T - 2*center->T + top->T)/dy/dy);
	//result->T = (left->T + 2*center->T + right->T + top->T + bottom->T)/6;

/*	// here should be navier-stokes equations

	// continuity equation (density)
	real dro = -dy*(right->density*right->vx - left->density*left->vx)/2/dx +
		(bottom->density*bottom->vy - top->density*top->vy)/2/dy;
	// momentum equation
	real dux_left = (center->vx - left->vx)/dx,
		dux_right = (right->vx - center->vx)/dx,
		duy_top = (center->vx - top->vx)/dx,
		duy_bottom = (bottom->vx - center->vx)/dx,
		duy_left
*/
	return 0;
}

real sqr(real x)
{
	return x*x;
}

real tau_xx(const real dux, const real dvy)
{
	return 4.f/3.f*dux - dvy/3.f;
}

real tau_xy(const real dux, const real duy, const real dvx, const real dvy)
{
	return duy + dvx - (dux + duy)/3.f;
}

real tau_yy(const real dux, const real dvy)
{
	return 4.f/3.f*dvy - dux/3.f;
}

int mech::solve_equation9(node *r,const edge9_t &edge,real dx,real dy,real dt) 
{
#if 0
	// simple thermal diffuse equation
	r->T = edge(0,0)->T + dt*((edge(1,0)->T - 2*edge(0,0)->T + edge(-1,0)->T)/dx/dx
		+ (edge(0,1)->T - 2*edge(0,0)->T + edge(0,-1)->T)/dy/dy);
#else
	const real Cp = R + Cv;
	// continuity equation (density)
	real rho = edge(0,0)->rho - dt*(
		(edge(1,0)->rho*edge(1,0)->u - edge(-1,0)->rho*edge(-1,0)->u)/dx +
		(edge(0,1)->rho*edge(0,1)->v - edge(0,-1)->rho*edge(0,-1)->v)/dy
		)/2.f;
	//real rho = edge(0,0)->rho;

	// momentum equation
	// various velocity derivatives
	real dux_left   = (edge( 0, 0)->u - edge(-1, 0)->u)/dx;
	real dux_right  = (edge( 1, 0)->u - edge( 0, 0)->u)/dx;
	real dux_top    = (edge( 1,-1)->u - edge(-1,-1)->u)/dx/2;
	real dux_bottom = (edge( 1, 1)->u - edge(-1, 1)->u)/dx/2;

	real duy_top    = (edge( 0, 0)->u - edge( 0,-1)->u)/dy;
	real duy_bottom = (edge( 0, 1)->u - edge( 0, 0)->u)/dy;
	real duy_left   = (edge(-1, 1)->u - edge(-1,-1)->u)/dy/2;
	real duy_right  = (edge( 1, 1)->u - edge( 1,-1)->u)/dy/2;

	real dvx_left   = (edge( 0, 0)->v - edge(-1, 0)->v)/dx;
	real dvx_right  = (edge( 1, 0)->v - edge( 0, 0)->v)/dx;
	real dvx_top    = (edge( 1,-1)->v - edge(-1,-1)->v)/dx/2;
	real dvx_bottom = (edge( 1, 1)->v - edge(-1, 1)->v)/dx/2;

	real dvy_top    = (edge( 0, 0)->v - edge( 0,-1)->v)/dy;
	real dvy_bottom = (edge( 0, 1)->v - edge( 0, 0)->v)/dy;
	real dvy_left   = (edge(-1, 1)->v - edge(-1,-1)->v)/dy/2;
	real dvy_right  = (edge( 1, 1)->v - edge( 1,-1)->v)/dy/2;

	// viscous stress
	real tau_xx_right = mu*(4.f/3.f*2.f*dux_right - 1.0f/3.0f*dvy_right);
	real tau_xx_left  = mu*(4.f/3.f*2.f*dux_left  - 1.0f/3.0f*dvy_left);

	real tau_xy_bottom = mu*(2.f*duy_bottom + dvx_bottom - 1.f/3.f*(dux_bottom + 2.f*dvy_bottom));
	real tau_xy_top    = mu*(2.f*duy_top    + dvx_top    - 1.f/3.f*(dux_bottom + 2.f*dvy_top));

	real tau_xy_right  = mu*(duy_right + 2*dvx_right - 1.f/3.f*(2*dux_right + dvy_right));
	real tau_xy_left   = mu*(duy_left  + 2*dvx_left  - 1.f/3.f*(2*dux_left  + dvy_left));

	//real tau_yx = mu*(2*(dvx_right - dvx_left) + duy_right - duy_left)
	//	- mu/3*(2*(dux_right - dux_left) + dvy_right - dvy_left);

	real tau_yy_bottom = mu*(4.f/3.f*2.f*dvy_bottom - 1.f/3.f*dux_bottom);
	real tau_yy_top    = mu*(4.f/3.f*2.f*dvy_top    - 1.f/3.f*dux_top);

	real momentum_udx = (
		+ edge(1,0)->rhou2() - edge(-1,0)->rhou2()
		+ edge(1,0)->p(R)    - edge(-1,0)->p(R)
		- tau_xx_right + tau_xx_left
		)/2/dx;

	real momentum_udy = (
		+ edge(0,1)->rhouv() - edge(0,-1)->rhouv()
		- tau_xy_bottom + tau_xy_top
		)/2/dy;

	real momentum_vdx = (
		+ edge(1,0)->rhouv() - edge(-1,0)->rhouv()
		- tau_xy_right + tau_xy_left
		)/2/dx;

	real momentum_vdy = (
		+ edge(0,1)->rhov2() - edge(0,-1)->rhov2()
		+ edge(0,1)->p(R)    - edge(0,-1)->p(R)
		- tau_yy_bottom + tau_yy_top
		)/2/dy;

	real rhou = edge(0,0)->rho*edge(0,0)->u - dt*(momentum_udx + momentum_udy); 
	real rhov = edge(0,0)->rho*edge(0,0)->v - dt*(momentum_vdx + momentum_vdy);

	// energy equation
	// diffuse
	real q_right  = (edge(1,0)->T - edge( 0,0)->T)/dx;
	real q_left   = (edge(0,0)->T - edge(-1,0)->T)/dx;
	real q_top    = (edge(0,0)->T - edge(0,-1)->T)/dy;
	real q_bottom = (edge(0,1)->T - edge( 0,0)->T)/dy;

	real energy_dx = (
		+ edge(1,0)->rho*edge(1,0)->u*edge(1,0)->e0(Cv) - edge(-1,0)->rho*edge(-1,0)->u*edge(-1,0)->e0(Cv)
		+ edge(1,0)->u*edge(1,0)->p(R) - edge(-1,0)->u*edge(-1,0)->p(R)
		- 2*lambda*(q_right - q_left)
		- edge(1,0)->u*tau_xx_right + edge(-1,0)->u*tau_xx_left
		- edge(1,0)->v*tau_xy_right + edge(-1,0)->v*tau_xy_left
		)/2/dx;

	real energy_dy = (
		+ edge(0,1)->rho*edge(0,1)->v*edge(0,1)->e0(Cv) - edge(0,-1)->rho*edge(0,-1)->v*edge(0,-1)->e0(Cv)
		+ edge(0,1)->v*edge(0,1)->p(R) - edge(0,-1)->v*edge(0,-1)->p(R)
		- 2*lambda*(q_bottom - q_top)
		- edge(0,1)->u*tau_xy_bottom + edge(0,-1)->u*tau_xy_top
		- edge(0,1)->v*tau_yy_bottom + edge(0,-1)->v*tau_yy_top
		)/2/dy;

	real e0 = edge(0,0)->e0(Cv) - dt*(energy_dx + energy_dy);

	if(rho < 0)rho = 0.0001f;
	r->rho = rho;
	r->u = rhou / rho;
	r->v = rhov / rho;
	r->T = (e0 - (r->u*r->u + r->v*r->v)/2)/Cv;
#endif
	if(r->T < 0 || r->rho < 0 || _isnan(r->T) || _isnan(r->rho) || _isnan(r->u) || _isnan(r->v)){
		//printf("bad value {%f,%f,%f,%f}\n",r->u,r->v,r->rho,r->T);
		//exit(-1);
		return -1;
	}
	return 0;
}