CS计算机代考程序代写 #include

#include
#include
#include
#include “datadef.h”
#include “init.h”

#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y)) extern int *ileft, *iright; extern int nprocs, proc; /* Computation of tentative velocity field (f, g) */ void compute_tentative_velocity(float **u, float **v, float **f, float **g, char **flag, int imax, int jmax, float del_t, float delx, float dely, float gamma, float Re) { int i, j; float du2dx, duvdy, duvdx, dv2dy, laplu, laplv; for (i=1; i<=imax-1; i++) { for (j=1; j<=jmax; j++) { /* only if both adjacent cells are fluid cells */ if ((flag[i][j] & C_F) && (flag[i+1][j] & C_F)) { du2dx = ((u[i][j]+u[i+1][j])*(u[i][j]+u[i+1][j])+ gamma*fabs(u[i][j]+u[i+1][j])*(u[i][j]-u[i+1][j])- (u[i-1][j]+u[i][j])*(u[i-1][j]+u[i][j])- gamma*fabs(u[i-1][j]+u[i][j])*(u[i-1][j]-u[i][j])) /(4.0*delx); duvdy = ((v[i][j]+v[i+1][j])*(u[i][j]+u[i][j+1])+ gamma*fabs(v[i][j]+v[i+1][j])*(u[i][j]-u[i][j+1])- (v[i][j-1]+v[i+1][j-1])*(u[i][j-1]+u[i][j])- gamma*fabs(v[i][j-1]+v[i+1][j-1])*(u[i][j-1]-u[i][j])) /(4.0*dely); laplu = (u[i+1][j]-2.0*u[i][j]+u[i-1][j])/delx/delx+ (u[i][j+1]-2.0*u[i][j]+u[i][j-1])/dely/dely; f[i][j] = u[i][j]+del_t*(laplu/Re-du2dx-duvdy); } else { f[i][j] = u[i][j]; } } } for (i=1; i<=imax; i++) { for (j=1; j<=jmax-1; j++) { /* only if both adjacent cells are fluid cells */ if ((flag[i][j] & C_F) && (flag[i][j+1] & C_F)) { duvdx = ((u[i][j]+u[i][j+1])*(v[i][j]+v[i+1][j])+ gamma*fabs(u[i][j]+u[i][j+1])*(v[i][j]-v[i+1][j])- (u[i-1][j]+u[i-1][j+1])*(v[i-1][j]+v[i][j])- gamma*fabs(u[i-1][j]+u[i-1][j+1])*(v[i-1][j]-v[i][j])) /(4.0*delx); dv2dy = ((v[i][j]+v[i][j+1])*(v[i][j]+v[i][j+1])+ gamma*fabs(v[i][j]+v[i][j+1])*(v[i][j]-v[i][j+1])- (v[i][j-1]+v[i][j])*(v[i][j-1]+v[i][j])- gamma*fabs(v[i][j-1]+v[i][j])*(v[i][j-1]-v[i][j])) /(4.0*dely); laplv = (v[i+1][j]-2.0*v[i][j]+v[i-1][j])/delx/delx+ (v[i][j+1]-2.0*v[i][j]+v[i][j-1])/dely/dely; g[i][j] = v[i][j]+del_t*(laplv/Re-duvdx-dv2dy); } else { g[i][j] = v[i][j]; } } } /* f & g at external boundaries */ for (j=1; j<=jmax; j++) { f[0][j] = u[0][j]; f[imax][j] = u[imax][j]; } for (i=1; i<=imax; i++) { g[i][0] = v[i][0]; g[i][jmax] = v[i][jmax]; } } /* Calculate the right hand side of the pressure equation */ void compute_rhs(float **f, float **g, float **rhs, char **flag, int imax, int jmax, float del_t, float delx, float dely) { int i, j; for (i=1;i<=imax;i++) { for (j=1;j<=jmax;j++) { if (flag[i][j] & C_F) { /* only for fluid and non-surface cells */ rhs[i][j] = ( (f[i][j]-f[i-1][j])/delx + (g[i][j]-g[i][j-1])/dely ) / del_t; } } } } /* Red/Black SOR to solve the poisson equation */ int poisson(float **p, float **rhs, char **flag, int imax, int jmax, float delx, float dely, float eps, int itermax, float omega, float *res, int ifull) { int i, j, iter; float add, beta_2, beta_mod; float p0 = 0.0; int rb; /* Red-black value. */ float rdx2 = 1.0/(delx*delx); float rdy2 = 1.0/(dely*dely); beta_2 = -omega/(2.0*(rdx2+rdy2)); /* Calculate sum of squares */ for (i = 1; i <= imax; i++) { for (j=1; j<=jmax; j++) { if (flag[i][j] & C_F) { p0 += p[i][j]*p[i][j]; } } } p0 = sqrt(p0/ifull); if (p0 < 0.0001) { p0 = 1.0; } /* Red/Black SOR-iteration */ for (iter = 0; iter < itermax; iter++) { for (rb = 0; rb <= 1; rb++) { for (i = 1; i <= imax; i++) { for (j = 1; j <= jmax; j++) { if ((i+j) % 2 != rb) { continue; } if (flag[i][j] == (C_F | B_NSEW)) { /* five point star for interior fluid cells */ p[i][j] = (1.-omega)*p[i][j] - beta_2*( (p[i+1][j]+p[i-1][j])*rdx2 + (p[i][j+1]+p[i][j-1])*rdy2 - rhs[i][j] ); } else if (flag[i][j] & C_F) { /* modified star near boundary */ beta_mod = -omega/((eps_E+eps_W)*rdx2+(eps_N+eps_S)*rdy2); p[i][j] = (1.-omega)*p[i][j] - beta_mod*( (eps_E*p[i+1][j]+eps_W*p[i-1][j])*rdx2 + (eps_N*p[i][j+1]+eps_S*p[i][j-1])*rdy2 - rhs[i][j] ); } } /* end of j */ } /* end of i */ } /* end of rb */ /* Partial computation of residual */ *res = 0.0; for (i = 1; i <= imax; i++) { for (j = 1; j <= jmax; j++) { if (flag[i][j] & C_F) { /* only fluid cells */ add = (eps_E*(p[i+1][j]-p[i][j]) - eps_W*(p[i][j]-p[i-1][j])) * rdx2 + (eps_N*(p[i][j+1]-p[i][j]) - eps_S*(p[i][j]-p[i][j-1])) * rdy2 - rhs[i][j]; *res += add*add; } } } *res = sqrt((*res)/ifull)/p0; /* convergence? */ if (*res= 1.0e-10) { /* else no time stepsize control */
umax = 1.0e-10;
vmax = 1.0e-10;
for (i=0; i<=imax+1; i++) { for (j=1; j<=jmax+1; j++) { umax = max(fabs(u[i][j]), umax); } } for (i=1; i<=imax+1; i++) { for (j=0; j<=jmax+1; j++) { vmax = max(fabs(v[i][j]), vmax); } } deltu = delx/umax; deltv = dely/vmax; deltRe = 1/(1/(delx*delx)+1/(dely*dely))*Re/2.0; if (deltu