Skip to content

Instantly share code, notes, and snippets.

@windstriver

windstriver/flow_over_flat_plate.c

Created December 23, 2014 12:52
Show Gist options
  • Save windstriver/e4e68ae01f67723916c9 to your computer and use it in GitHub Desktop.
Save windstriver/e4e68ae01f67723916c9 to your computer and use it in GitHub Desktop.
Download ZIP
N-S equation, time marching
Raw
flow_over_flat_plate.c
/*
* Navierstokes Time marching 2D.
* Balaji Sankar.
*/
//The Only stuff you need to Change are in this box.
//*Please*Do not tamper with the other stuff unless you know exactly what you are doing.
//Blasius solution
//=============================================================
// Number of divisions in X axis
#define IMAX 70
// Number of divisions in Y axis
#define JMAX 90
// Number of iterations to be done.
#define CALLMAC 1000
// Number of iterations after which intermediate results are printed.
#define MACPRINT 1000
// Enter 1 to get a very long output file with all the nitty gritties displayed.
#define VERBOSITY 0
//==============================================================
#include<stdio.h>
#include<math.h>
// Defining geometrical constants. 10 cm long plate,1mm separation
#define LHORI .1
#define LVERT (.001*LHORI)
#define R 287
#define gamma 1.4
#define cv (R/(gamma-1))
#define cp (R*gamma/(gamma-1))
#define t_inf 288
#define p_inf 1.01325e5
//Your wish
#define Mach_inf 0.2
#define Mach_wall 0.01
#define t_wall 300
#define a_inf (pow ( ( gamma*R*t_inf ),0.5 ))
//Your wish
#define u_inf (a_inf*Mach_inf)
//Your wish
#define v_inf (a_inf*0)
//Density
#define rho_inf (p_inf/(R*t_inf))
// Reynolds number based on plate length
#define Re_plate (rho_inf*u_inf*LHORI/VISC(t_inf))
//Comes often.
#define IJ (IMAX*JMAX)
//Prandtl's number. I have no idea why it is 0.71
#define Pr 0.71
//define modes for computing txx ,tyy and txy
// prediction in x direction , while in e
#define XPRED 1
// prediction in y direction , while in f
#define YPRED 2
// correction in x direction , while in e
#define XCORR 3
// correction in y direction , while in f
#define YCORR 4
// Macro calculate viscosity using Sutherlands formula,reference conditions are 1.789e-5 at 288.16 k.
#define VISC(t) (1.789e-5* ( pow ( ( t/288.16 ), ( 3/2 ) ) ) * ( ( 288.16+110 ) / ( t+110 ) ))
//Grid spacing. Dont change manually. Change IMAX,JMAX if required.
#define dx (LHORI/ ( IMAX-1 ))
#define dy (LVERT/ ( JMAX-1 ))
//Function prototype declaration so that we can return non int values to calling function.
void bc( float*);
void mac(float*,FILE*,FILE*,float*,float* ,float* ,float* ,float* ,float* ,float* ,float* ,float*,int);
void fp2u(float*,float*);
void fp2e(float*,float*,int);
void fp2f(float*,float*,int);
float txx(float*,int,int);
float tyy(float*,int,int);
float txy(float*,int,int);
float qx(float* ,int,int);
float qy(float* ,int,int);
void dudt_predictor(float* ,float* ,float*,float* );
void dudt_corrector(float* ,float* ,float*,float* );
void u_predict(float * ,float * ,float *,float);
void u_correct(float * ,float * ,float * ,float *,float);
void utofp(float * ,float *);
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
int main()
{
printf("dx=%e dy=%e ",dx,dy);
// Creating flow properties matrixmac(float * afp,FILE *fid1,float * au,float * ae,float * af,float * aup,float * adudt,float * adudtp,float * afpp)
float fp[5][IMAX][JMAX]={0},u[5][IMAX][JMAX]={0},e[5][IMAX][JMAX]={0},f[5][IMAX][JMAX]={0},ep[5][IMAX][JMAX]={0},fp_vector[5][IMAX][JMAX]={0},up[5][IMAX][JMAX]={0},dudt[5][IMAX][JMAX]={0},dudtp[5][IMAX][JMAX]={0},fpp[5][IMAX][JMAX]={0};
// iteration variables stored in cpu register for speed
register int i,j,k;
// Opening file to dump output while debugging
FILE *fid1,*fid2;
//fid1=fopen("output.txt","w");
fid1=fopen("output.csv","w");
fid2=fopen("u_profile.txt","w");
// Assign initial conditions in interior points only
for ( k=0;k<=4;k++ )
{
for ( i=1;i<= ( IMAX-2 );i++ )
{
for ( j=1;j<= ( JMAX-2 );j++ )
{
if (k==0) fp[k][i][j]=u_inf;
if (k==1) fp[k][i][j]=v_inf;
if (k==2) fp[k][i][j]=p_inf;
if (k==3) fp[k][i][j]=t_inf;
if (k==4) fp[k][i][j]=fp[k-2][i][j]/ ( 287*fp[k-1][i][j] );
}
}
}
//===========================================================
// Call bc assigner to get fp on borders;
bc(&fp[0][0][0]);
//===========================================================
//Print fp b4 callin MAC
fprintf(fid1,"\n Printing Initial flow properties in main, Boundaries assigned. \n");
for ( k=0;k<=4;k++ )
{
fprintf(fid1,"Flow property number %d \n \n",k+1);
for ( i=0;i<= ( IMAX-1 );i++ )
{
for ( j=0;j<= ( JMAX-1 );j++ )
{
fprintf(fid1," %e, ",fp[k][i][j]);
}
fprintf(fid1," \n ");
}
fprintf(fid1," \n ");
}
//===========================================================
// Call mac how many times?
if (VERBOSITY==1)// call verbose version
{
for (i=1;i<=CALLMAC;i++)
verbose_mac(&fp[0][0][0],fid1,u,e,f,ep,fp_vector,up,dudt,dudtp,fpp,i);
}
else
{
for (i=1;i<=CALLMAC;i++)
mac(&fp[0][0][0],fid1,fid2,u,e,f,ep,fp_vector,up,dudt,dudtp,fpp,i);
}
//===========================================================
//Print fp after callin MAC
fprintf(fid1,"\n Printing Values from main, Output. \n Govinda Govinda Govindaaaaaa........ \n\n");
for ( k=0;k<=4;k++ )
{
switch (k)
{
case(0):
fprintf(fid1,"Velocity Parallel to plate.Flow property number %d \n",k+1);
break;
case(1):
fprintf(fid1,"Velocity Perpendicular to plate.Flow property number %d \n",k+1);
break;
case(2):
fprintf(fid1,"Static Pressure .Flow property number %d \n",k+1);
break;
case(3):
fprintf(fid1,"Static Temperature .Flow property number %d \n",k+1);
break;
case(4):
fprintf(fid1,"Density .Flow property number %d \n",k+1);
break;
}
for ( i=0;i<= ( IMAX-1 );i++ )
{
for ( j=0;j<= ( JMAX-1 );j++ )
{
fprintf(fid1," %e, ",fp[k][i][j]);
}
fprintf(fid1," \n ");
}
fprintf(fid1," \n ");
}
//===========================================================
printf ("\nGee.. Reload the output file doc :-) \n\n\n\n");
//===========================================================
//writing a u file for gnu plot.
/*for ( j=1;j<=JMAX;j++ )
{
fprintf(fid2,"%d %f \n ",j,fp[0][IMAX-1][j-1]);
}
*/fclose(fid2);
fclose(fid1);
return(0);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining bc function in ANSI c style with type declaration in parameter reception itself.
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void bc(float *afp)
{
register int i,j;
//===========================================================
//Leading Edge
*(afp+(0)*IJ)=0; //u
*(afp+(1)*IJ)=0; //v, next page , first locationfor leadin edge.
*(afp+(2)*IJ)=p_inf; //p
*(afp+(3)*IJ)=t_inf; //t
*(afp+(4)*IJ)=p_inf/(R*t_inf); //density
//===========================================================
//INflow,j starts from 1 as LE is assigned above.
//Top row of each page; that is each property's top row is assigned.left side is wall.
for ( j=1;j<=JMAX-1;j++ )
{
*(afp+(0)*IJ+j) =u_inf;
*(afp+(1)*IJ+j) =v_inf;
*(afp+(2)*IJ+j) =p_inf;
*(afp+(3)*IJ+j) =t_inf;
*(afp+(4)*IJ+j) =p_inf/(R*t_inf);
}
//===========================================================
//Upper surface, Right side of matrix
for ( i=2;i<=IMAX;i++ )
{
*(afp+(0)*IJ+(i*JMAX)-1) =u_inf;
*(afp+(1)*IJ+(i*JMAX)-1) =v_inf;
*(afp+(2)*IJ+(i*JMAX)-1) =p_inf;
*(afp+(3)*IJ+(i*JMAX)-1) =t_inf;
*(afp+(4)*IJ+(i*JMAX)-1) =p_inf/(R*t_inf);
}
//===========================================================
//Outflow surface
for ( j=1;j<JMAX-1;j++ )
{
*(afp+(0)*IJ+j+(IMAX-1)*(JMAX)) =(2* *(afp+(0)*IJ+j+(IMAX-2)*(JMAX)))-*(afp+(0)*IJ+j+(IMAX-3)*(JMAX));
*(afp+(1)*IJ+j+(IMAX-1)*(JMAX)) =(2* *(afp+(1)*IJ+j+(IMAX-2)*(JMAX)))-*(afp+(1)*IJ+j+(IMAX-3)*(JMAX));
*(afp+(2)*IJ+j+(IMAX-1)*(JMAX)) =p_inf;
*(afp+(3)*IJ+j+(IMAX-1)*(JMAX)) =(2* *(afp+(3)*IJ+j+(IMAX-2)*(JMAX)))-*(afp+(3)*IJ+j+(IMAX-3)*(JMAX));
*(afp+(4)*IJ+j+(IMAX-1)*(JMAX)) = *(afp+(2)*IJ+j+(IMAX-1)*(JMAX))/(R* *(afp+(3)*IJ+j+(IMAX-1)*(JMAX)));
}
//===========================================================
//Base Plate;
for ( i=1;i<=IMAX-1;i++ )
{
*(afp+(0)*IJ+i*JMAX) =0;
*(afp+(1)*IJ+i*JMAX) =0;
*(afp+(2)*IJ+i*JMAX) =( 2* *(afp+(2)*IJ+i*JMAX+1))-(*(afp+(2)*IJ+i*JMAX+2));
*(afp+(3)*IJ+i*JMAX) =t_wall;
*(afp+(4)*IJ+i*JMAX) =*(afp+(2)*IJ+i*JMAX)/(R*t_wall);
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining mac function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void mac(float * afp,FILE *fid1,FILE *fid2,float * au,float * ae,float * af,float * aep,float * afp_vector,float * aup,float * adudt,float * adudtp,float * afpp,int maccall)
{
register int i,j,k;
float v_dash_max=-100000,v_dash,dt_min=1000,dt,redx,redy,redxmax=-1000000,redxmin=1000000,redymax=-1000000,redymin=100000;
//Simplifying expression for dt
float k1=(1/(dx*dx))+(1/(dy*dy));
float a,ubydx,vbydy;
for (i=1;i<=IJ;i++)
{
v_dash=(4*(pow(VISC(*(afp+i-1+3*IJ)),2)*gamma)/(Pr*3* *(afp+i-1+4*IJ)));
redx=*(afp+i-1+4*IJ) * *(afp+i-1+0*IJ) *dx/VISC(*(afp+i-1+3*IJ));
redx=*(afp+i-1+4*IJ) * *(afp+i-1+1*IJ) *dy/VISC(*(afp+i-1+3*IJ));
//Finding out the max value of v_dash
if (v_dash>=v_dash_max)
{
v_dash_max=v_dash;
}
if (redx>=redxmax)
{
redxmax=redx;
}
if (redy>=redymax)
{
redymax=redy;
}
if (redx<=redxmin)
{
redxmin=redx;
}
if (redy<=redymin)
{
redymin=redy;
}
}
// fprintf(fid1,"\n");
for (i=1;i<=IJ;i++)
{
//Sound velocity
a=(sqrt(gamma*R* *(afp+i-1+3*IJ)));
ubydx=(fabs(*(afp+(i-1)+(0*IJ)))/dx);
vbydy=(fabs(*(afp+(i-1)+(1*IJ)))/dy);
dt=1/(ubydx+vbydy+(a*sqrt(k1))+(2*v_dash_max*k1));
if (dt<=dt_min)
{
dt_min=dt;
}
}
dt=dt_min*0.8;
fprintf(fid1,"\ndx=,%e ,dy=,%e ,dt =, %e, i=, %d",dx,dy,dt_min,maccall);
fprintf(fid1,"\nredxmax=,%e ,redxmin=,%e ,redymax =, %e, redymin=, %e",redxmax,redxmin,redymax,redymin);
//Populate u
fp2u(au,afp);
//populate dudt (flux vector) in prediction mode.
dudt_predictor(adudt,ae,af,afp);
//populate U_predicted in prediction mode.
u_predict(adudt,au,aup,dt);
//populate fp from up_predict
u2fp(aup,afpp);
//Call to bc to apply bc to fpp
bc(afpp);
//Compute the future dudt using predicted flow properties, e anf f will be changed in the process.
dudt_corrector(adudtp,aep,afp_vector,afpp);
//compute corrected u , replacing the old u.
u_correct(adudt,adudtp,au,aup,dt);
//populate fp from up_corrected
u2fp(au,afp);
bc(afp);
//Printing current values every 100 iterations or so , for safety.
if (maccall%MACPRINT==0)
{
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Velocity Parallel to plate \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+0*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Velocity Perpendicular to plate \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+1*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Pressure \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+2*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Temperature \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+3*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Density \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+4*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Velocity Parallel to plate \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+0*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Velocity Perpendicular to plate \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+1*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Pressure \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+2*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Temperature \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+3*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing Intermediate results from MAC at iteration number= %d. Density \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+4*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
//Writing Output file for blasius flow
for ( j=0;j<=JMAX-1;j++ )
{
fprintf(fid2,"%f %d \n ",*(afp+0*IJ+((IMAX-1)*JMAX)+j),j+1);
}
}
}
void verbose_mac(float * afp,FILE *fid1,float * au,float * ae,float * af,float * aep,float * afp_vector,float * aup,float * adudt,float * adudtp,float * afpp,int maccall)
{
register int i,j,k;
float v_dash_max=-10,v_dash,dt_min=100,dt,redx,redy,redxmax=-1000,redxmin=1000,redymax=-10000,redymin=1000;
//Simplifying expression for dt
float k1=(1/(dx*dx))+(1/(dy*dy));
float a,ubydx,vbydy;
for (i=1;i<=IJ;i++)
{
v_dash=(4*(pow(VISC(*(afp+i-1+3*IJ)),2)*gamma)/(Pr*3* *(afp+i-1+4*IJ)));
//Finding out the max value of v_dash
if (v_dash>=v_dash_max)
{
v_dash_max=v_dash;
}
if (redx>=redxmax)
{
redxmax=redx;
}
if (redy>=redymax)
{
redymax=redy;
}
if (redx<=redxmin)
{
redxmin=redx;
}
if (redy<=redymin)
{
redymin=redy;
}
}
// fprintf(fid1,"\n");
for (i=1;i<=IJ;i++)
{
//Sound velocity
a=(sqrt(gamma*R* *(afp+i-1+3*IJ)));
ubydx=(fabs(*(afp+(i-1)+(0*IJ)))/dx);
vbydy=(fabs(*(afp+(i-1)+(1*IJ)))/dy);
dt=1/(ubydx+vbydy+(a*sqrt(k1))+(2*v_dash_max*k1));
if (dt<=dt_min)
{
dt_min=dt;
}
}
dt=dt_min*0.8;
fprintf(fid1,"\ndx=,%e ,dy=,%e ,dt =, %e, i=, %d",dx,dy,dt_min,maccall);
//fprintf(fid1,"\nredxmax=,%e ,redxmin=,%e ,redymax =, %e, redymin=, %e",redxmax,redxmin,redymax,redymin);
//display flow properties as received by MAC
fprintf(fid1,"\n Printing as received fp at iteration number= %d. Velocity Parallel to plate \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+0*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing as received fp at iteration number= %d. Velocity Perpendicular to plate \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+1*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing as received fp at iteration number= %d. Pressure \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+2*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing as received fp at iteration number= %d. Temperature \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+3*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
fprintf(fid1,"\n Printing as received fp at iteration number= %d. Density \n",maccall);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",(*(afp+i-1+4*IJ)));
if ((i%(JMAX)==0)&&(i!=0))
{
//Printing a new line once JMAX is reached.
fprintf(fid1," \n");
}
}
//Populate u
fp2u(au,afp);
//populate dudt (flux vector) in prediction mode.
dudt_predictor(adudt,ae,af,afp);
fprintf(fid1,"\n Printing Values from MAC. \n Predicted values of dudt FROM CURRENT Flow Properties... \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing predicted dudt number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(adudt+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//printing e and f
fprintf(fid1,"\n Printing Values from MAC. \n e used in dudt computation... \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing e from current number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(ae+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
fprintf(fid1,"\n Printing Values from MAC. \n f used in dudt computation... \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing f from current number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(af+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//populate U_predicted in prediction mode.
u_predict(adudt,au,aup,dt);
fprintf(fid1,"\n Printing Values from MAC. \nPredicted values of U FROM dudt FROM CURRENT Floe Properties... \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing Predicted values of U number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(aup+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//populate fp from up_predict
u2fp(aup,afpp);
fprintf(fid1,"\n Printing Values from MAC. \nPredicted values of fp,fpp AT THE NEXT TIME STEP... \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing Predicted flow property number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(afpp+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//Call to bc to apply bc to fpp
bc(afpp);
//Compute the future dudt using predicted flow properties, e anf f will be changed in the process.
fprintf(fid1,"\n Printing Values from MAC. \nPredicted values of fp,fpp. after call 2 bc.. \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing Predicted values of floe properties after BC, number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(afpp+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
dudt_corrector(adudtp,aep,afp_vector,afpp);
fprintf(fid1,"\n Printing Values from MAC. \n dudt using predicted fp \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing dudt computed using predicted fp, number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(adudtp+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//compute corrected u , replacing the old u.
u_correct(adudt,adudtp,au,aup,dt);
fprintf(fid1,"\n Printing Values from MAC. \n Corrected U... \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing Corrected U using averaged dudt, number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(au+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//populate fp from up_corrected
u2fp(au,afp);
bc(afp);
fprintf(fid1,"\n Printing Values from MAC. \nfp from corrected u \n\n");
for (k=0;k<=3;k++)
{
fprintf(fid1,"\n Printing Final flow property number %d \n",k+1);
for (i=1;i<=IJ;i++)
{
fprintf(fid1," %e,",*(afp+i-1+k*IJ));
if ((i%(JMAX)==0)&&(i!=0))
{
fprintf(fid1," \n");
}
}
}
//Printing current values every 100 iterations or so , for safety.
if (maccall%MACPRINT==0)
{
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining fp2u function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void fp2u(float * au,float * afp)
{
int i;
for (i=1;i<=IJ;i++)
{
*(au+i-1+0*IJ)=*(afp+i-1+4*IJ); //density
*(au+i-1+1*IJ)=*(afp+i-1+4*IJ) * *(afp+i-1+0*IJ); //density * u
*(au+i-1+2*IJ)=*(afp+i-1+4*IJ) * *(afp+i-1+1*IJ); //density * v
*(au+i-1+3*IJ)=*(afp+i-1+4*IJ) * ((cv* *(afp+i-1+3*IJ))+(pow(*(afp+i-1+0*IJ),2)+(pow(*(afp+i-1+1*IJ),2)))/2); //energy rho*(cv t+(u^2+v^2)/2);
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining fp2e function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void fp2e(float * ae,float * afp,int mode)
{
int i;
for (i=1;i<=IJ;i++)
{
*(ae+i-1+0*IJ)=*(afp+i-1+4*IJ) * *(afp+i-1+0*IJ) ; //density* u
*(ae+i-1+1*IJ)=(*(afp+i-1+4*IJ) * pow(*(afp+i-1+0*IJ),2))+(*(afp+i-1+2*IJ))-(txx(afp,i,mode)); //rho*u^2 +p -txx
*(ae+i-1+2*IJ)=(*(afp+i-1+4*IJ) * *(afp+i-1+0*IJ) * *(afp+i-1+1*IJ))-(txy(afp,i,mode)); //rho*u*v -txy
*(ae+i-1+3*IJ)=((*(afp+i-1+4*IJ)*((cv* *(afp+i-1+3*IJ))+((pow(*(afp+i-1+0*IJ),2)+pow(*(afp+i-1+1*IJ),2))/2)))+*(afp+i-1+2*IJ)) * *(afp+i-1+0*IJ)-(*(afp+i-1+0*IJ)*(txx(afp,i,mode)))-(*(afp+i-1+1*IJ)*(txy(afp,i,mode)))+(qx(afp,i,mode)); // u(p+rho*e)-u*txx-v*txy+qx
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining fp2f function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void fp2f(float * af,float * afp,int mode)
{
int i;
for (i=1;i<=IJ;i++)
{
*(af+i-1+0*IJ)=*(afp+i-1+4*IJ) * *(afp+i-1+1*IJ) ; //density* v
*(af+i-1+1*IJ)=(*(afp+i-1+4*IJ) * *(afp+i-1+0*IJ) * *(afp+i-1+1*IJ))-(txy(afp,i,mode)); //rho*u*v -txy
*(af+i-1+2*IJ)=(*(afp+i-1+4*IJ) * pow(*(afp+i-1+1*IJ),2))+(*(afp+i-1+2*IJ))-(tyy(afp,i,mode)); //rho*v^2 +p -tyy
float et=(*(afp+i-1+4*IJ)*((cv* *(afp+i-1+3*IJ))+((pow(*(afp+i-1+0*IJ),2)+pow(*(afp+i-1+1*IJ),2))/2)));
float p=*(afp+i-1+2*IJ);
float v=*(afp+i-1+1*IJ);
float u=*(afp+i-1+0*IJ);
*(af+i-1+3*IJ)=((et+p)*v)-(v*tyy(afp,i,mode))-(u*txy(afp,i,mode))+qy(afp,i,mode); // u(p+rho*e)-v*tyy-u*txy+qy
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining txx function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
float txx(float *afp,int i,int mode)
{
float mu =VISC(*(afp+i-1+3*IJ));
float lambda =(mu)*(-2/3);
float txx_r, dvdy,dudx;
if (mode==XPRED)//prediction in x dir'
{//perpendicular(Y) direction has central difference and along(X) uses in rear ward difference in prediction mode
if ((i-1)%JMAX==0)
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ+1) - *(afp+i-1+1*IJ));//On flat plate forward difference is the only option
}
else if ((i)%JMAX==0)
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ) - *(afp+i-1+1*IJ-1));//On top surface backward difference is the only possibility
}
else
{
dvdy=(0.5/dy)*(*(afp+i-1+1*IJ+1) - *(afp+i-1+1*IJ-1));//On the interior central difference is used
}
if (i<=JMAX)
{
dudx=(1/dx)*(*(afp+i-1+0*IJ+JMAX)-*(afp+i-1+0*IJ));//On inlet forward difference
}
else
{
dudx=(1/dx)*(*(afp+i-1+0*IJ)-*(afp+i-1+0*IJ-JMAX));//%Every where else (including exit) rear ward difference.
}
}
else if (mode==XCORR)//correction in x dir
{//perpendicular(Y) direction has central difference and along(X) uses in forward difference in correction mode
if ((i-1)%JMAX==0)
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ+1) - *(afp+i-1+1*IJ));//On flat plate forward difference is the only option
}
else if ((i)%JMAX==0)
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ) - *(afp+i-1+1*IJ-1));//On top surface backward difference is the only possibility
}
else
{
dvdy=(0.5/dy)*(*(afp+i-1+1*IJ+1) - *(afp+i-1+1*IJ-1));//On the interior central difference is used
}
if (i>(JMAX*(IMAX-1)))
{
dudx=(1/dx)*(*(afp+i-1+0*IJ )-*(afp+i-1+0*IJ-JMAX));; //On exit rearward difference
}
else
{
dudx=(1/dx)*(*(afp+i-1+0*IJ+JMAX)-*(afp+i-1+0*IJ));//Every where else (including inlet) forward difference.
}
}
else
{
printf("\n\n Some problem in mode declaration in txx \n");
}
txx_r=lambda*(dudx+dvdy)+2*mu*(dudx);
return(txx_r);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining tyy function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
float tyy(float *afp,int i,int mode)
{
float mu =VISC(*(afp+i-1+3*IJ));
float lambda =(mu)*(-2/3);
float tyy_r, dvdy,dudx;
if (mode==YPRED)//prediction in y dir'
{//perpendicular(X) direction has central difference and along(Y) uses rear ward difference in prediction mode
if ((i-1)%JMAX==0)
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ+1) - *(afp+i-1+1*IJ));//On flat plate forward difference is the only option
}
else
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ) - *(afp+i-1+1*IJ-1));//On interior rearward difference in y direction
}
if (i<=JMAX)
{
dudx=(1/dx)*(*(afp+i-1+0*IJ+JMAX)-*(afp+i-1+0*IJ));//On inlet forward difference
}
else if (i>(JMAX*(IMAX-1)))
{
dudx=(1/dx)*(*(afp+i-1+0*IJ)-*(afp+i-1+0*IJ-JMAX));//%Exit rear ward difference.
}
else
{
dudx=(0.5/dx)*(*(afp+i-1+0*IJ+JMAX) - *(afp+i-1+0*IJ-JMAX));//On the interior central difference is used
}
}
else if (mode==YCORR)//correction in y dir
{//perpendicular(x) direction has central difference and along(y) uses in forward difference in correction mode
if (i<=JMAX)
{
dudx=(1/dx)*(*(afp+i-1+0*IJ+JMAX)-*(afp+i-1+0*IJ));//On inlet forward difference
}
else if (i>(JMAX*(IMAX-1)))
{
dudx=(1/dx)*(*(afp+i-1+0*IJ)-*(afp+i-1+0*IJ-JMAX));//%Exit rear ward difference.
}
else
{
dudx=(0.5/dx)*(*(afp+i-1+0*IJ+JMAX) - *(afp+i-1+0*IJ-JMAX));//On the interior central difference is used
}
if ((i)%JMAX==0)
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ) - *(afp+i-1+1*IJ-1));//On top rear difference is the only option
}
else
{
dvdy=(1/dy)*(*(afp+i-1+1*IJ+1) - *(afp+i-1+1*IJ));//On interior forward difference in y direction
}
}
else
{
printf("\n\n Some problem in mode declaration in tyy \n");
}
tyy_r=lambda*(dudx+dvdy)+2*mu*(dvdy);
return(tyy_r);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining txy function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
float txy(float *afp,int i,int mode)
{
float mu =VISC(*(afp+i-1+3*IJ));
float lambda =(mu)*(-2/3);
float txy_r, dvdx,dudy;
if (mode==XPRED)//prediction in x dir'
{//perpendicular(Y) direction has central difference and along(X) uses in rear ward difference in prediction mode
if ((i-1)%JMAX==0)
{
dudy=(1/dy)*(*(afp+i-1+0*IJ+1) - *(afp+i-1+0*IJ));//On flat plate forward difference is the only option
}
else if ((i)%JMAX==0)
{
dudy=(1/dy)*(*(afp+i-1+0*IJ) - *(afp+i-1+0*IJ-1));//On top surface backward difference is the only possibility
}
else
{
dudy=(0.5/dy)*(*(afp+i-1+0*IJ+1) - *(afp+i-1+0*IJ-1));//On the interior central difference is used
}
if (i<=JMAX)
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ+JMAX)-*(afp+i-1+1*IJ));//On inlet forward difference
}
else
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ)-*(afp+i-1+1*IJ-JMAX));//%Every where else (including exit) rear ward difference.
}
}
else if (mode==XCORR) //correction in x direction
{//perpendicular(Y) direction has central difference and along(X) uses in forward difference in correction mode
if ((i-1)%JMAX==0)
{
dudy=(1/dy)*(*(afp+i-1+0*IJ+1) - *(afp+i-1+0*IJ));//On flat plate forward difference is the only option
}
else if ((i)%JMAX==0)
{
dudy=(1/dy)*(*(afp+i-1+0*IJ) - *(afp+i-1+0*IJ-1));//On top surface backward difference is the only possibility
}
else
{
dudy=(0.5/dy)*(*(afp+i-1+0*IJ+1) - *(afp+i-1+0*IJ-1));//On the interior central difference is used
}
if (i>(JMAX*(IMAX-1)))
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ )-*(afp+i-1+1*IJ-JMAX));; //On exit rearward difference
}
else
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ+JMAX)-*(afp+i-1+1*IJ));//Every where else (including inlet) forward difference.
}
}
else if (mode==YPRED)//prediction in y dir'
{//perpendicular(X) direction has central difference and along(Y) uses rear ward difference in prediction mode
if ((i-1)%JMAX==0)
{
dudy=(1/dy)*(*(afp+i-1+0*IJ+1) - *(afp+i-1+0*IJ));//On flat plate forward difference is the only option
}
else
{
dudy=(1/dy)*(*(afp+i-1+0*IJ) - *(afp+i-1+0*IJ-1));//On interior rearward difference in y direction
}
if (i<=JMAX)
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ+JMAX)-*(afp+i-1+1*IJ));//On inlet forward difference
}
else if (i>(JMAX*(IMAX-1)))
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ)-*(afp+i-1+1*IJ-JMAX));//%Exit rear ward difference.
}
else
{
dvdx=(0.5/dx)*(*(afp+i-1+1*IJ+JMAX) - *(afp+i-1+1*IJ-JMAX));//On the interior central difference is used
}
}
else if (mode==YCORR)//correction in y dir
{//perpendicular(x) direction has central difference and along(y) uses in forward difference in correction mode
if (i<=JMAX)
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ+JMAX)-*(afp+i-1+1*IJ));//On inlet forward difference
}
else if (i>(JMAX*(IMAX-1)))
{
dvdx=(1/dx)*(*(afp+i-1+1*IJ)-*(afp+i-1+1*IJ-JMAX));//%Exit rear ward difference.
}
else
{
dvdx=(0.5/dx)*(*(afp+i-1+1*IJ+JMAX) - *(afp+i-1+1*IJ-JMAX));//On the interior central difference is used
}
if ((i)%JMAX==0)
{
dudy=(1/dy)*(*(afp+i-1+0*IJ) - *(afp+i-1+0*IJ-1));//On top rear difference is the only option
}
else
{
dudy=(1/dy)*(*(afp+i-1+0*IJ+1) - *(afp+i-1+0*IJ));//On interior forward difference in y direction
}
}
else
{
printf("\n\n Some problem in mode declaration in txy \n");
}
txy_r=mu*(dudy+dvdx);
return(txy_r);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining qx function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
float qx(float *afp,int i,int mode )
{
float mu=VISC(*(afp+i-1+3*IJ));
float k=mu*cp/Pr;
float dtdx,qx_r;
if (mode==1)
{
if (i<=JMAX)
{
dtdx=(1/dx)*(*(afp+i-1+3*IJ+JMAX)-*(afp+i-1+3*IJ));//On inlet forward difference
}
else
{
dtdx=(1/dx)*(*(afp+i-1+3*IJ)-*(afp+i-1+3*IJ-JMAX));//%Every where else (including exit) rear ward difference.
}
}
else if (mode==3)
{
if (i>(JMAX*(IMAX-1)))
{
dtdx=(1/dx)*(*(afp+i-1+3*IJ )-*(afp+i-1+3*IJ-JMAX)); //On exit rearward difference
}
else
{
dtdx=(1/dx)*(*(afp+i-1+3*IJ+JMAX)-*(afp+i-1+3*IJ)); //Every where else (including inlet) forward difference.
}
}
else
printf("Some problem in mode declaration in qx \n");
qx_r=-k*dtdx;
return(qx_r);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining qy function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
float qy(float *afp,int i,int mode )
{
float mu=VISC(*(afp+i-1+3*IJ));
float k=mu*cp/Pr;
float dtdy,qy_r;
// printf("\n k=%e",k);
if (mode==2)
{//perpendicular(X) direction has central difference and along(Y) uses rear ward difference in prediction mode
if ((i-1)%JMAX==0)
{
dtdy=(1/dy)*(*(afp+i-1+3*IJ+1) - *(afp+i-1+3*IJ));//On flat plate forward difference is the only option
}
else
{
dtdy=(1/dy)*(*(afp+i-1+3*IJ) - *(afp+i-1+3*IJ-1));//On interior rearward difference in y direction
}
}
else if (mode==4)//correction in y dir
{//perpendicular(x) direction has central difference and along(y) uses in forward difference in correction mode
if ((i)%JMAX==0)
{
dtdy=(1/dy)*(*(afp+i-1+3*IJ) - *(afp+i-1+3*IJ-1));//On top rear difference is the only option
}
else
{
dtdy=(1/dy)*(*(afp+i-1+3*IJ+1) - *(afp+i-1+3*IJ));//On interior forward difference in y direction
}
}
qy_r=-k*dtdy;
return(qy_r);
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining dudt_predict function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void dudt_predictor(float * dudt,float * ae,float *af,float *afp)
{
int i;
//use prediction mode to populate e,f vectors
fp2e(ae,afp,1);//derivative in x direction, prediction
fp2f(af,afp,2);//derivative in y direction, prediction
for (i=1;i<=IJ;i++)
{
if ( (i<=JMAX) || ((i-1)%JMAX==0) ||((i)%JMAX==0) || ((i>(JMAX*(IMAX-1)))) )
{
continue;
}
*(dudt+i-1+0*IJ)=-(1/dx)*(*(ae+i-1+0*IJ+JMAX)-*(ae+i-1+0*IJ))-(1/dy)*(*(af+i-1+0*IJ+1) - *(af+i-1+0*IJ));
*(dudt+i-1+1*IJ)=-(1/dx)*(*(ae+i-1+1*IJ+JMAX)-*(ae+i-1+1*IJ))-(1/dy)*(*(af+i-1+1*IJ+1) - *(af+i-1+1*IJ));
*(dudt+i-1+2*IJ)=-(1/dx)*(*(ae+i-1+2*IJ+JMAX)-*(ae+i-1+2*IJ))-(1/dy)*(*(af+i-1+2*IJ+1) - *(af+i-1+2*IJ));
*(dudt+i-1+3*IJ)=-(1/dx)*(*(ae+i-1+3*IJ+JMAX)-*(ae+i-1+3*IJ))-(1/dy)*(*(af+i-1+3*IJ+1) - *(af+i-1+3*IJ));
*(dudt+i-1+4*IJ)=0;
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining dudt_correct function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void dudt_corrector(float * dudtp,float * aep,float *afp_vector,float *afpp)
{
int i;
//use prediction mode to populate e,f vectors
fp2e(aep,afpp,3);//derivative in x direction, correction
fp2f(afp_vector,afpp,4);//derivative in y direction, correction
for (i=1;i<=IJ;i++)
{
if ( (i<=JMAX) || ((i-1)%JMAX==0) ||((i)%JMAX==0) || ((i>(JMAX*(IMAX-1)))) )
{
continue;
}
*(dudtp+i-1+0*IJ)=-(1/dx)*(*(aep+i-1+0*IJ)-*(aep+i-1+0*IJ-JMAX))-(1/dy)*(*(afp_vector+i-1+0*IJ) - *(afp_vector+i-1+0*IJ-1));
*(dudtp+i-1+1*IJ)=-(1/dx)*(*(aep+i-1+1*IJ)-*(aep+i-1+1*IJ-JMAX))-(1/dy)*(*(afp_vector+i-1+1*IJ) - *(afp_vector+i-1+1*IJ-1));
*(dudtp+i-1+2*IJ)=-(1/dx)*(*(aep+i-1+2*IJ)-*(aep+i-1+2*IJ-JMAX))-(1/dy)*(*(afp_vector+i-1+2*IJ) - *(afp_vector+i-1+2*IJ-1));
*(dudtp+i-1+3*IJ)=-(1/dx)*(*(aep+i-1+3*IJ)-*(aep+i-1+3*IJ-JMAX))-(1/dy)*(*(afp_vector+i-1+3*IJ) - *(afp_vector+i-1+3*IJ-1));
*(dudtp+i-1+4*IJ)=0;
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining u_predict function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void u_predict(float *adudt,float *au,float *aup,float dt)
{
int i;
for (i=1;i<=IJ;i++)
{
if ( (i<=JMAX) || ((i-1)%JMAX==0) ||((i)%JMAX==0) || ((i>(JMAX*(IMAX-1)))) )
{
continue;
}
*(aup+i-1+0*IJ)= dt* (*(adudt+i-1+0*IJ)) + *(au+i-1+0*IJ);
*(aup+i-1+1*IJ)= dt* (*(adudt+i-1+1*IJ)) + *(au+i-1+1*IJ);
*(aup+i-1+2*IJ)= dt* (*(adudt+i-1+2*IJ)) + *(au+i-1+2*IJ);
*(aup+i-1+3*IJ)= dt* (*(adudt+i-1+3*IJ)) + *(au+i-1+3*IJ);
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining u_correct function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void u_correct(float*adudt,float*adudtp,float*au,float*aup,float dt)
{
int i;
for (i=1;i<=IJ;i++)
{
if ( (i<=JMAX) || ((i-1)%JMAX==0) ||(i%JMAX==0) || ((i>(JMAX*(IMAX-1)))) )
{
continue;
}
float temp1= 0.5*(*(au+i-1+0*IJ)+*(aup+i-1+0*IJ)+dt*(*(adudtp+i-1+0*IJ)));
float temp2= 0.5*(*(au+i-1+1*IJ)+*(aup+i-1+1*IJ)+dt*(*(adudtp+i-1+1*IJ)));
float temp3= 0.5*(*(au+i-1+2*IJ)+*(aup+i-1+2*IJ)+dt*(*(adudtp+i-1+2*IJ)));
float temp4= 0.5*(*(au+i-1+3*IJ)+*(aup+i-1+3*IJ)+dt*(*(adudtp+i-1+3*IJ)));
*(au+i-1+0*IJ)=temp1;
*(au+i-1+1*IJ)=temp2;
*(au+i-1+2*IJ)=temp3;
*(au+i-1+3*IJ)=temp4;
}
}
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//Defining u2fp function in ANSI c style,
//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
void u2fp(float* aup,float*afpp)
{
int i;
for (i=1;i<=IJ;i++)
{
if ( (i<=JMAX) || ((i-1)%JMAX==0) ||((i)%JMAX==0) || ((i>(JMAX*(IMAX-1)))) )
{
continue;
}
*(afpp+0*IJ+i-1)=*(aup+i-1+1*IJ)/ *(aup+i-1+0*IJ);//U_predicted(i,j,2)/U_predicted(i,j,1);
*(afpp+1*IJ+i-1)=*(aup+i-1+2*IJ)/ *(aup+i-1+0*IJ);//U_predicted(i,j,3)/U_predicted(i,j,1);
*(afpp+3*IJ+i-1)=((*(aup+i-1+3*IJ)/ *(aup+i-1+0*IJ))-(pow(*(afpp+i-1+1*IJ),2)+pow(*(afpp+i-1+0*IJ),2))/2)/cv;//;((U_predicted(i,j,4)/U_predicted(i,j,1))-((fp_predicted(i,j,1))^2+(fp_predicted(i,j,2))^2)/2)/cv;
*(afpp+2*IJ+i-1)=(*(aup+i-1+0*IJ)*R* *(afpp+i-1+3*IJ));
*(afpp+4*IJ+i-1)=*(aup+i-1+0*IJ);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment