/*********************************************************************/
/* Title : Toyama CoAstal flow Simulator (TCAS) */
/*********************************************************************/
/* */
/* Author: Shoichi FURUYAMA */
/* */
/* Affiliation : Toyama National College of Maritime Technology */
/* JAPAN */
/* */
/* - Visiting Scholar at University of Queensland */
/* Australia (March 2007 - March 2008) */
/* */
/* E-mail: shoichi@toyama-cmt.ac.jp */
/* */
/* WWW: http://www.toyama-cmt.ac.jp/~shoichi/ */
/* */
/* Version: Ver. TCAS20080104A */
/* */
/* Development Environment: */
/* - OS: Redhat Linux 2.4.21-32.EL */
/* - Compiler: Intel(R) C Compiler for */
/* Intel(R) EM64T-based applications, Version 8.1 */
/* */
/* Copyright 2008, Shoichi FURUYAMA */
/*********************************************************************/
/*********************************************************************
.......................................................................
..... Conditions of the calculation ...................................
.......................................................................
0.30m +-------Open boundary --------+
(24)  |                              |
      |                              |
      |                           <- air flow
0.0875m +-----------------------------+ 0.0875m(7)
       |                          <- inlet
0.0125m|                    Water <- inlet
(1)      <- outlet                <- inlet
 1.89    +-----------------------------+
   [m/s]                           5.625m(450)
+ dx=dy=0.0125 [m]
+ IMAX=450 JMAX=24
+ Smagorinsky Coefficient: Cs=0.1
+ MUwater 0.001002           water viscosity [Pa s]
+ RHOwater 998.2             density of water [kg/m^3]
+ RHOair 1.166               density of air [kg/m^3]
+ dt=5.0e-4 [s]
........................................................................
..... Visualization (examples for visualization on GNUPLOT) ............
........................................................................
.. file ..
’./Data/avs????.dat’
.. Free-surface profile (example) ..
gnuplot> set size ratio -1
gnuplot> set dgrid3d 50,80
gnuplot> set pm3d map
gnuplot> splot "avs0002.dat" every::2 using 1:2:7
.. vector profile (example) ..
gnuplot> vs=0.1
gnuplot> plot "avs0004.dat" every::2 using 1:2:($5*vs):($6*vs) w vec
**********************************************************************/
#include<stdio.h>
#include<math.h>
#include<stdlib.h>//slx.add
#define DX 0.0125
#define DY 0.0125
#define IMAX 450
#define JMAX 28
/* INLET & OUTLET geometry. The number of meshes, NOT LENGTH! */
#define INLET 7
#define OUTLET 1
#define TMAX 20. /* [sec] */
#define STPMAX 5 /* in the STPMAX=1000 case of dt=5.0e-4, the end of calculation is T=0.5 */
#define RT 0.1
#define OUTPUT 40 /* In the case of 50[Hz], OUTPUT*dt is set to 0.02(=2.0*10e-2).
If dt is 10e-5, the similar condition is that the OUTPUT is set to 2000,
If dt is 10e-4, the similar condition is that the OUTPUT is set to 200,
=> If dt is 5.0*10e-4, the similar condition is that the OUTPUT is set to 40,
If dt is 10e-3, the similar condition is that the OUTPUT is set to 20. */
#define RCOUTP (OUTPUT*10)
#define OUTPUT_ITE 2
#define OUTPUT_CAV 200
#define OUTP_MIN 2000 /* Change OUTPUT to 1 */
#define ALPHA 0.5
#define BETA 1.0
#define EPS (1.e-12) /* residual for SOR */
#define C99 0.99
#define MITE 10000 /* maximum iteration count for SOR */
#define GAMMA 1.4 /* specific heat rate */
#define KAPPA 7.15 /* experimental value */
#define B (3049.*100000) /* [Pa] experimental value */
#define MUair (1.82e-5) /* air viscosity [Pa s] */
#define MUwater (1.002e-3) /* water viscosity [Pa s] */
#define VKC 0.4 /* von Kerman const. */
#define Cs 0.1 /* Smagorinsky const. 0.1 - 0.2*/
#define CVIS 0.0 /* artificial viscosity effection coef.*/
#define Patm 101325. /* pressure of atmosphere [Pa] or [N/m^2] or [kg.m/s^2] */
#define RHOair 1.166 /* density of air [kg/m^3] */
#define RHOwater 998.2 /* density of water [kg/m^3] */
#define G 9.8 /* gravitational acc. */
#define TANFAC 0.9
#define RNOISE 0.8
#define CP 50
#define PI 3.141592
/* INLET condition */
#define IN_RHO RHOwater
#define IN_PHI 1.
#define IN_U -1.0
#define V0 IN_U
#define D0 (INLET*DY)
#define IN_V 0.
#define OUT_U -1.89
#define OUT_V 0.
#define N 7 /* smooth turbulence */
#define DD0R 0.75 /* delta / d0 */
#define Vmax (V0/( 1.-( DD0R/(N+1) ) ) )
/* monitorling line x */
#define X1 17
#define X2 167
#define X3 213
#define X4 263
#define X5 313
#define X6 413
double rho[IMAX][JMAX], gxrho[IMAX][JMAX], gzrho[IMAX][JMAX], rhon[IMAX][JMAX], gxrhon[IMAX][JMAX], gzrhon[IMAX][JMAX];
double u[IMAX][JMAX], gxu[IMAX][JMAX], gzu[IMAX][JMAX], un[IMAX][JMAX], gxun[IMAX][JMAX], gzun[IMAX][JMAX];
double w[IMAX][JMAX], gxw[IMAX][JMAX], gzw[IMAX][JMAX], wn[IMAX][JMAX], gxwn[IMAX][JMAX], gzwn[IMAX][JMAX];
double phi[IMAX][JMAX], gxphi[IMAX][JMAX], gzphi[IMAX][JMAX], phin[IMAX][JMAX], gxphin[IMAX][JMAX], gzphin[IMAX][JMAX];
double p[IMAX][JMAX], gxp[IMAX][JMAX], gzp[IMAX][JMAX], pn[IMAX][JMAX], gxpn[IMAX][JMAX], gzpn[IMAX][JMAX];
double u_s[IMAX][JMAX],w_s[IMAX][JMAX],u_w[IMAX][JMAX],w_u[IMAX][JMAX];



enum true_false{TRUE=-2, FALSE};
double T, DT, ITVLM;
double max(double x, double y){
return( (x > y) ? x : y);
}
double min(double x, double y){
return( (x < y) ? x : y);
}
double set_dt(double u[IMAX][JMAX], double v[IMAX][JMAX]){
double dt;
dt=5.0e-4;
return(dt);
}
double dphi_digit(double phi){
//double dphi(double phi){
// phi (0<= phi <=1 ) is transformed to tangental space value dphi */
double ddphi;
// if(phi<0. || phi>1.){
// printf("(dphi)INVALID phi value(=%16.8e)\n",phi);
// exit(-1);
// }
ddphi=tan(C99*PI*(phi-0.5));
return(ddphi);
}
//double dphi_ND(double phi){
double dphi(double phi){
/* for nondigitizing method, same value returning.. */
// phi (0<= phi <=1 ) is transformed to tangental space value dphi */
double ddphi;
// if(phi<0. || phi>1.){
	// printf("(dphi)INVALID phi value(=%16.8e)\n",phi);
// exit(-1);
// }
// ddphi=tan(C99*PI*(phi-0.5));
ddphi=phi;
return(ddphi);
}
double iphi_digit(double phi){
//double iphi(double phi){
/* digitizing */
double iiphi;
// phi (tangental space) is transformed to normal value (0<= iphi <=1) */
iiphi=atan(phi)/(C99*PI)+0.5;
// if(iiphi>1. || iiphi<0.){
// printf("(iphi) INVALID iiphi value(=%16.8e)\n",iiphi);
// exit(-1);
// }
return(iiphi);
}
double iphi(double phi){
//double iphi_ND(double phi){
/* for non-digitizing, only return same value */
double iiphi;
iiphi=phi;
return(iiphi);
}
int velo_calc(double u[IMAX][JMAX], double v[IMAX][JMAX],
double u_s[IMAX][JMAX], double v_s[IMAX][JMAX],
double u_v[IMAX][JMAX], double v_u[IMAX][JMAX]){
int i,j;
/*
u_s : u at scalar point
v_s : v at scalar point
u_v : u at velocity v point
v_u : v at velocity u point
*/
for(i=0; i<=IMAX-2; i++){
for(j=0; j<=JMAX-2; j++){
u_s[i][j]=0.5*(u[i][j]+u[i+1][j]);
v_s[i][j]=0.5*(v[i][j]+v[i][j+1]);
}
}
for(i=1; i<=IMAX-2; i++){
for(j=1; j<=JMAX-2; j++){
u_v[i][j]=0.25*(u[i][j]+u[i+1][j ]+u[i ][j-1]+u[i+1][j-1]);
v_u[i][j]=0.25*(v[i][j]+v[i ][j+1]+v[i-1][j ]+v[i-1][j+1]);
}
}
return(0);
}
int advection_rcip(double f[IMAX][JMAX], double fx[IMAX][JMAX], double fy[IMAX][JMAX],
double u[IMAX][JMAX], double v[IMAX][JMAX], double fn[IMAX][JMAX]){
/* rcip.... from sample */
int i, j, isx, jsy;
double a1, a2, a3, a4, a5, a6, a7, a8, a9;
double b1, b2, b3, b4;
double c1, c2, cx1, cy2;
double xi, yi, fxj, fyi;
double au, av;
double ddx, ddy, fxdd, fydd, rxd, ryd;
double fxn[IMAX][JMAX], fyn[IMAX][JMAX];
double sx, sy, theta;
double oldf;
theta=1.;
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
fn[i][j]=f[i][j];
fxn[i][j]=fx[i][j];
fyn[i][j]=fy[i][j];
}
}
for(i=1; i<IMAX-1; i++){
for(j=1; j<JMAX-1; j++){
au=u[i][j];
av=v[i][j];
rxd=1.;
ryd=1.;
if(au <= 0.0){
ddx=DX;
isx=1;
}
else{
ddx=-DX;
isx=-1;
}
if(av <= 0.0){
ddy=DY;
jsy=1;
}
else{
ddy=-DY;
jsy=-1;
}
xi=-au*DT;
yi=-av*DT;
sx=(fn[i+isx][j]-fn[i][j])/ddx;
fxj=fxn[i+isx][j]-sx;
if(fabs(fxj)<=1.0e-10){
fxj=1.0e-10;
fxn[i][j]=sx;
}
sy=(fn[i][j+jsy]-fn[i][j])/ddy;
fyi=fyn[i][j+jsy]-sy;
if(fabs(fyi)<=1.0e-10){
fyi=1.0e-10;
fyn[i][j]=sy;
}
fxdd=fxn[i+isx][j]*fxn[i][j];
fydd=fyn[i][j+jsy]*fyn[i][j];
b1=sx-fxn[i][j];
b2=fxj;
c1=(fabs(b1/b2)-1.)/ddx;
if(fxdd <= 0.0) c1=c1*theta*rxd;
else c1=0.0;
cx1=1.+c1*ddx;
b3=sy-fyn[i][j];
b4=fyi;
c2=(fabs(b3/b4)-1.)/ddy;
if(fydd <= 0.0) c2=c2*theta*ryd;
else c2=0.0;
cy2=1.+c2*ddy;
a1=fxn[i][j]+c1*fn[i][j];
a2=fyn[i][j]+c2*fn[i][j];
a7=(cy2*fyn[i][j+jsy]-sy*cy2+fyn[i][j]-sy)/ddy/ddy;
a6=(cx1*fxn[i+isx][j]-sx*cx1+fxn[i][j]-sx)/ddx/ddx;
a9=(cy2*fn[i][j+jsy]-fn[i][j]-a2*ddy)/ddy/ddy-a7*ddy;
a8=(cx1*fn[i+isx][j]-fn[i][j]-a1*ddx)/ddx/ddx-a6*ddx;
a3=cy2*fxn[i][j+jsy]/ddy+cx1*fyn[i+isx][j]/ddx
+c1*fn[i][j+jsy]/ddy+c2*fn[i+isx][j]/ddx
+(fn[i][j]-(1.+c1*ddx+c2*ddy)*fn[i+isx][j+jsy])/ddx/ddy
+a6*ddx*ddx/ddy+a7*ddy*ddy/ddx+a8*ddx/ddy+a9*ddy/ddx;
a5=cx1*fyn[i+isx][j]/ddx/ddx+c2*fn[i+isx][j]/ddx/ddx
-a2/ddx/ddx-a3/ddx;
a4=cy2*fxn[i][j+jsy]/ddy/ddy+c1*fn[i][j+jsy]/ddy/ddy
-a1/ddy/ddy-a3/ddy;
f[i][j]=fn[i][j]+a1*xi+a2*yi+a3*xi*yi+a4*xi*yi*yi
+a5*xi*xi*yi+a6*xi*xi*xi+a7*yi*yi*yi+a9*yi*yi
+a8*xi*xi;
if((1.+c1*xi+c2*yi)==0.){
printf("0 divide in advection c1=%8.4f c2=%8.4f xi=%8.4f yi=%8.4f at (%d,%d)\n",c1,c2,xi,yi,i,j);
exit(-1);
}
f[i][j]=f[i][j]/(1.+c1*xi+c2*yi);
fx[i][j]=a1+a3*yi+a4*yi*yi+2.*a5*xi*yi+3.*a6*xi*xi
+2.*a8*xi-c1*f[i][j];
fx[i][j]=fx[i][j]/(1.+c1*xi+c2*yi);
fy[i][j]=a2+a3*xi+a5*xi*xi+2.*a4*xi*yi+3.*a7*yi*yi
+2.*a9*yi-c2*f[i][j];
fy[i][j]=fy[i][j]/(1.+c1*xi+c2*yi);
}
}
for(i=2; i<IMAX-1; i++){
for(j=2; j<JMAX-1; j++){
oldf=fn[i][j];
fn[i][j]=f[i][j];
f[i][j]=oldf;
fxn[i][j]=fx[i][j]
-fx[i][j]*(u[i+1][j]-u[i-1][j])*0.5*DT/DX
-fy[i][j]*(v[i+1][j]-v[i-1][j])*0.5*DT/DX;
fyn[i][j]=fy[i][j]
-fx[i][j]*(u[i][j+1]-u[i][j-1])*0.5*DT/DY
-fy[i][j]*(v[i][j+1]-v[i][j-1])*0.5*DT/DY;
}
}
return(0);
}
int advection(double f[IMAX][JMAX], double GX[IMAX][JMAX], double GY[IMAX][JMAX],
double u[IMAX][JMAX], double v[IMAX][JMAX], double fn[IMAX][JMAX]){
/* normal CIP */
int i, j, isgn, jsgn, im1, jm1;
double XX, YY;
double A1, A2, A3, A4, A5, A6, A7, A8, TMP;
double GXt[IMAX][JMAX], GYt[IMAX][JMAX];
for(i=1; i<IMAX-1; i++){
for(j=1; j<JMAX-1; j++){
XX=-u[i][j]*DT;
YY=-v[i][j]*DT;
if(u[i][j]>=0.)isgn=1;
else isgn=-1;
if(v[i][j]>=0)jsgn=1;
else jsgn=-1;
im1=i-isgn;
jm1=j-jsgn;
A8=f[i][j]-f[im1][j]-f[i][jm1]+f[im1][jm1];
TMP=GY[im1][j]-GY[i][j];
A1=((GX[im1][j]+GX[i][j])*DX*isgn-2.*(f[i][j]-f[im1][j]))/(DX*DX*DX*isgn);
A2=(-A8-(GX[i][jm1]-GX[i][j])*DX * isgn)/(DX*DX*DY*jsgn);
A3=(3.*(f[im1][j]-f[i][j])+(GX[im1][j]+2.*GX[i][j])*DX*isgn)/(DX*DX);
A4=(A2*DX*DX-TMP)/(DX*isgn);
A5=(-2*(f[i][j]-f[i][jm1])+(GY[i][jm1]+GY[i][j])*DY*jsgn)/(DY*DY*DY*jsgn);
A6=(-A8-TMP*DY*jsgn)/(DX*DY*DY*isgn);
A7=(3.*(f[i][jm1]-f[i][j])+(GY[i][jm1]+2.*GY[i][j])*DY*jsgn)/(DY*DY);
fn[i][j]=((A1*XX+A2*YY+A3)*XX+A4*YY+GX[i][j])*XX
+((A5*YY+A6*XX+A7)*YY+GY[i][j])*YY+f[i][j];
GX[i][j]=(3.*A1*XX+2.*(A2*YY+A3))*XX + (A4+A6*YY)*YY+GX[i][j];
GY[i][j]=(3.*A5*YY+2.*(A6*XX+A7))*YY + (A4+A2*XX)*XX+GY[i][j];
}
}
return(0);
}
int err_output(double p[IMAX][JMAX], double rho[IMAX][JMAX], double u[IMAX][JMAX], double v[IMAX][JMAX], double phi[IMAX][JMAX],
int count, double time){
// for microAVS
FILE *fpfld, *fpdat;
int i, j;
int stp;
char avsfld[25]= "./Data/err????.fld";
char avsdat[25]= "./Data/err????.dat";
char avsfldl[25]= "./err????.fld";
char avsdatl[25]= "./err????.dat";
stp=count/OUTPUT;
avsfld[10] = '0'   + stp/1000;
avsfld[11] = '0' + (stp%1000)/100;
avsfld[12] = '0' + ((stp%1000)%100)/10;
avsfld[13] = '0' + ((stp%1000)%100)%10;
avsdat[10] = '0' + stp/1000;
avsdat[11] = '0' + (stp%1000)/100;
avsdat[12] = '0' + ((stp%1000)%100)/10;
avsdat[13] = '0' + ((stp%1000)%100)%10;
avsdatl[5] = '0' + stp/1000;
avsdatl[6] = '0' + (stp%1000)/100;
avsdatl[7] = '0' + ((stp%1000)%100)/10;
avsdatl[8] = '0' + ((stp%1000)%100)%10;
avsfldl[5] = '0' + stp/1000;
avsfldl[6] = '0' + (stp%1000)/100;
avsfldl[7] = '0' + ((stp%1000)%100)/10;
avsfldl[8] = '0' + ((stp%1000)%100)%10;
if ((fpfld=fopen (avsfld,"w+")) == NULL){
printf("%s was not opened!\n",avsfld);
printf("’Maybe you should make ’Data’ directory....\n");
exit(1);
}
fprintf(fpfld,"# AVS field file\n");
fprintf(fpfld,"#\n");
fprintf(fpfld,"ndim = 2\n");
fprintf(fpfld,"dim1 = %d\n",JMAX);
fprintf(fpfld,"dim2 = %d\n",IMAX);
fprintf(fpfld,"nspace = 2\n");
fprintf(fpfld,"veclen = 5\n");
fprintf(fpfld,"data = double\n");
fprintf(fpfld,"field = irregular\n");
fprintf(fpfld,"label = pressure density u v phi\n");
fprintf(fpfld,"\n");
fprintf(fpfld,"coord 1 file=%s filetype=ascii skip=2 offset=0 stride=7\n",avsdatl);
fprintf(fpfld,"coord 2 file=%s filetype=ascii skip=2 offset=1 stride=7\n",avsdatl);
fprintf(fpfld,"variable 1 file=%s filetype=ascii skip=2 offset=2 stride=7\n",avsdatl);
fprintf(fpfld,"variable 2 file=%s filetype=ascii skip=2 offset=3 stride=7\n",avsdatl);
fprintf(fpfld,"variable 3 file=%s filetype=ascii skip=2 offset=4 stride=7\n",avsdatl);
fprintf(fpfld,"variable 4 file=%s filetype=ascii skip=2 offset=5 stride=7\n",avsdatl);
fprintf(fpfld,"variable 5 file=%s filetype=ascii skip=2 offset=6 stride=7\n",avsdatl);
fclose(fpfld);
if ((fpdat=fopen (avsdat,"w+")) == NULL){
printf("%s was not opened!\n",avsdat);
printf("’Maybe you should make ’Data’ directory....\n");
exit(1);
}
fprintf(fpdat,"t=%16.8e\n",time);
fprintf(fpdat,"%16s %16s %16s %16s %16s %16s %16s\n","x","y","pressure","density","u","w","phi");
for(i=0; i<IMAX; i++){
for(j=0; j<JMAX; j++){
fprintf(fpdat,"%16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
(i-2)*DX,(j-2)*DY,p[i][j],rho[i][j],u[i][j],v[i][j],iphi(phi[i][j]));
}
}
fclose(fpdat);
return(0);
}
int output1(double f[IMAX][JMAX], int count){
FILE *fp;
int i, j;
int stp;